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Q '. Abstract 

We present a simple formalism for the evolution of time like jets in which tree-level matrix 
element corrections can be systematically incorporated, up to arbitrary parton multiplicities and 
i over all of phase space, in a way that exponentiates the matching corrections. The scheme is 

cast as a shower Markov chain which generates one single unweighted event sample, that can 
P> |- be passed to standard hadronization models. Remaining perturbative uncertainties are estimated 

^ ' by providing several alternative weight sets for the same events, at a relatively modest additional 

overhead. As an explicit example, we consider Z qq evolution with unpolarized, massless 
quarks and include several formally subleading improvements as well as matching to tree-level 
matrix elements through a^. The resulting algorithm is implemented in the publicly available 
ViNCiA plugin' to the PythiaS event generator. 
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Csl '. 1 Introduction 

(N 

O . The experimental program now underway at the Large Hadion Collider (LHC) will make extensive 

demands on theorists' ability to predict background processes, governed by Standard-Model physics. 
These predictions require perturbative computations in the electroweak sector, and call on both pertur- 
bative and nonperturbative physics in quantum chromodynamics. The dependence on nonfactorized 
^ I and nonperturbative QCD is not yet amenable to a first-principles calculation, and must therefore rely 

on experimental measurements of the parton distribution functions, as well as models for hadroniza- 
tion and the underlying event. An appropiiate choice of observables — infrared- and coUinear-safe 
ones — can minimize (but not eliminate) the dependence on these models. The factorizable perturba- 
tive component of the backgrounds can be computed systematically. 

An essential class of observables for new-physics searches are multi-jet differential cross sec- 
tions, typically with events also required to include decay products of one or more electroweak vector 
bosons. Jet shapes and jet substructure observables are also important, both from a calibration point 
of view as well as to search for decays of boosted objects [1,2]. The two basic approaches to com- 
puting perturbative contributions are via a fixed-order expansion in powers of the strong coupling Og, 
and in a parton-shower approach which resums leading (LL) and possibly next-to-leading logarithms 
(NLL) of large ratios of scales multiplying the strong coupling. The former approach can be carried 
out using widely-available tools to leading order (LO) in ag for essentially any jet multiplicity, for 
a growing list of processes to next-to-leading order (NLO), and for a select list to next-to-next-to- 
leading order (NNLO). It sacrifices a detailed picture of each event, and of jet substructure, in favor of 
of a systematically improved description of the "hard" or wide-angle radiation reflected in the several 
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jets. The parton-shower approach (see [3] for a recent review) favors developing an accurate picture 
of the "soft" or "collinear" radiation that dresses the hard event, thereby providing an event-by-event 
description of jet substmcture as well as allowing the incorporation of a nonperturbative model for the 
final-state hadrons. 

From a theoretical point of view, the two approaches correspond to summing up different but 
overlapping sets of contributions to the short-distance perturbative matrix elements. Simply adding 
the two would yield an overcounting of the common contributions, and the separation of the two 
sets is somewhat delicate. The recent decade has seen the appearance of a number of strategies 
for combining or "matching" these two approaches for general final states. These can be broadly 
described as "slicing", "subtraction", or "unitarity" approaches. 

In a slicing approach, the phase space for multiple emissions is separated into two regions. In 
one, the calculation uses the leading-order matrix element; in the other, the leading-log approxima- 
tion as given by a parton shower. The schemes introduced by Mangano (MLM) in the days of the top 
quark discovery and fomialized later [4,5], and by Catani, Krauss, Kuhn, and Webber (CKKW) [6] 
ai^e examples of a slicing approach. The CKKW approach has been implemented within the Sherpa 
framework [7]; the MLM one, using Alpgen [8] interfaced [9, 10] to both Pythia [11] and Her- 
WIG [12]. Refined versions of the CKKW method have since been introduced by Lonnblad [13, 14] 
in the context of the color-dipole model [15, 16] and implemented in the Ariadne generator [17], 
as well as by Mrenna and Richardson [18] using MadGraph [19], again interfaced to Herwig and 
Pythia. The original strategy [20, 21] for matching in Herwig, for one emission beyond a basic 
process, also follows this approach and may be seen as a precursor to the CKKW formalism. 

In a subtractive approach, the shower approximation is subtracted from the exact fixed-order ma- 
trix element, and two classes of events are generated: standard events and counter-events. The latter 
may have negative weights. The Mc@NLO program built on Herwig is an example of subtrac- 
tive matching to one loop [22]. Other subtraction implementations include a Sherpa implementa- 
tion [23] and one by Dinsdale, Ternick, and Weinzierl (DTW) [24] based on the Catani-Seymour 
(CS) method [25] of fixed-order calculations. Finally, the Menlops scheme [26] probably repre- 
sents the most advanced current matching method, combining POWHEG (a unitarity-based variant of 
Mc@NLO, see below) with a slicing -based matching for multijet emissions. 

Note, however, that in both slicing and subtraction, any subleading divergencies in the matched 
matrix elements beyond one additional parton are not regulated by the (LL) shower, and hence all the 
multileg schemes — MLM, CKKW, and MENLOPS — are forced to introduce a "matching scale" 
below which only the pure LL shower is used. While such schemes therefore guarantee the rates of 
hard additional jets to be correct to LO, the same is not guaranteed for jet substructure, which can be 
explicitly sensitive to multiparton correlations below the matching scale and beyond LL [27]. 

In a unitarity approach, to maintain a sequence of unit-weight events, the selection of branching 
events is modified by a veto depending on the ratio of the exact matrix element to the shower approxi- 
mation. Since the correction is applied on the splitting probability itself, it is automatically resummed 
to all orders by the shower Sudakovs, with real and virtual corrections canceling order by order in per- 
turbation theory. Indeed, the original approach to matching, caixied out for one additional emission 
beyond a basic process in Pythia [28,29], follows this approach. It is also used in Powheg [30], 
there combined with a subtractive matching to NLO. The proposal of matching by Sudakov reweight- 
ing by Nagy and Soper [31] is also within this class. However, the unitary approach to matching has 
so far only been worked out for a single emission. We shall here generalize it to an arbitraiy number 
of emissions, aixiving at the equivalent to Menlops but based on unitarity instead of slicing for the 
additional emissions. This will allow us to extend the matching over all of phase space, and should 
therefore result in a more accurate modeling not only of jet rates but also of jet substructure. 
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Most showers, with the exception of Ariadne and the Winter-Krauss shower [32], are based on 
collinear factorization, which is to say 1 — 2 branching in shower evolution. (Pythia 8 combines 
a 1 — )• 2 splitting probability with a 2 — )• 3 phase-space mapping.) In the present paper, we continue 
the development of a leading-log (LL) parton shower [33] based on dipole antennae, that is 2 — 3 
branching. We choose a simpler context than hadron collisions, that of electron-positron collisions. 
This allows us to set aside the questions of initial-state emission as well as those of the underlying 
event. 

In sec. 2, we describe in greater detail the ingredients needed for such a shower, as well as our 
normalization conventions, and compai^e the origins of different singularities and corresponding log- 
arithms in different shower formalisms. We also discuss the different matching approaches in more 
detail. In sec. 3, we discuss the evolution integral, and show how to cast it in a general form whose 
specializations correspond to a wide variety of interesting evolution variables. We then solve the re- 
sulting evolution equation. In sec. 4, we discuss the shower algorithm, as well as improvements that 
can be made to its logarithmic accuracy. In sec. 5, we discuss the details of matching the dipole- 
antenna shower to tree-level matrix elements, at both leading and subleading color. The procedure 
we use to evaluate the remaining perturbative uncertainties is described in sec. 6, and in sec. 7, we 
comment on hadronization; in sec. 8, we compare the results of running the unitarity-based approach 
implemented in ViNCiA to LEP data and to PYTHIA 8. We make some concluding remarks in sec. 9. 



2 Nomenclature and Conventions 

In this section, we introduce the basic elements of our perturbative formalism, which is largely based 
on ref. [33]. First, in sec. 2.1, we illustrate how the KLN theorem may be used to rewrite the coeffi- 
cients of perturbation theory as the expansion of an all-orders Markov chain, using NLO as an explicit 
example. Then, in sec. 2.2, we briefly describe each of the ingredients that enter our dipole-antenna 
shower formalism. 



2.1 Perturbation Theory with Markov Chains 

Consider the Bom-level cross section for an arbitrary hard process, H, differentially in an arbitrary 
infrared-safe observable O, 



dan 



dO 



Bom 



Jd^H \mP\^6{0-0{{p}h)), (1) 



where the integration mns over the full final-state on-shell phase space of H (this expression and 
those below would also apply to hadron collisions were we to include integrations over the parton 
distribution functions in the initial state), and the 5 function projects out a 1-dimensional slice defined 
by O evaluated on the set of final-state momenta which we denote {p}h (without the 5 function, the 
integration over phase space would just give the total cross section, not the differential one). 

To make the connection to parton showers, and to discuss all-orders resummations in that context, 
we may insert an operator, S, that acts on the Born-level final state before the observable is evaluated, 
i.e., 

dan 



dO 



S 



Jd<!>H \mP\^S{{p}h,0). (2) 



Formally, this operator — the evolution operator — will be responsible for generating all (real and 
virtual) higher-order conections to the Born-level expression. The measurement 6 function appear- 
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ing explicitly in eq. (1) is now implicit in S. (Ultimately, non-perturbative corrections can also be 
included.) 

Algorithmically, we shall cast S as an iterative Markov chain, with an evolution parameter that 
formally represents the factorization scale of the event, below which all structure is summed over 
inclusively. As the Markov chain develops, the evolution parameter will go towards zero, and the 
event structure will become more and more exclusively resolved. A transition from a perturbative 
evolution to a non-perturbative one can also be inserted, at an appropriate scale, typically around 
1 GeV. This scale thus represents the lowest perturbative scale that can appear in the calculations, 
with all perturbative corrections below it summed over inclusively. 

It is instinctive to begin by considering the first-order expansion the operator must have in order 
to agree with NLO perturbation theory. 



with the one-loop amplitude and the ratio d^u+i /^^^H in the second Une representing the 
phase space of one additional final-state paiticle; we shall return to the associated factorization be- 
low. The two coiTcction teniis ai^e separately divergent and hence eq. (3) only has a symbolic formal 
meaning. It requires a regulator for actual calculations. Introducing the factorization scale mentioned 
above, and introducing an n + 1 — )• n mapping of momenta by summing inclusively over all emissions 
below it, we obtain, instead, the first-order expansion corresponding to an evolution from the starting 
scale, s (the cm. energy squai^ed), down to the scale Q\, 

where the factorization scale, Qe (a.k.a. the "evolution scale"), separates resolved from unresolved 
regions. This expression is well-defined if the functional form of Q^; properly separates singular 
from non-singular regions, i.e., is "infrared sensible" [34]. (Coixections to this expression arising 
from scales below Qe will be taken into account by eventually letting Qe — )• 0.) Due to the KLN 
theorem [35], the real and virtual singularities must be equal and of opposite sign, thus we can rewrite 

|M(?)p Jo d<^H |m(?¥ ' 

where is a non-singular function when the regulator is removed, allowing us to express eq. (4) 
as 

S^'H{p}H,s,Ql,0) = (i + kP- r^^^-^%lf]6iO-0{{p}H)) 

^4 dC |,g.^ '°-°(M«.0). (« 
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In this form, the NLO correction to the total cross section is given solely by the term Kjj:\ with the 
remaining terms having been written in an explicitly unitai^y construction. 

We can also rewrite the exact ratio |M//+ip/|M//p as a process-dependent term whose integral 
is non-singulai; plus a sum over universal singular- ones, 

where r runs over "radiators", whose precise definition, such as partons or dipoles, depends on the 
chosen decomposition of the singular- structures in jAfjif+ip, and the superscript [''1 on the phase 
space factors indicate that each radiator may in principle be associated with a different phase space 
factorization. 

By the simple rewritings above, we have now obtained a form of the expansion in which the singu- 
larity and unitarity stiiicture of S are both explicitly manifest. Deviations from unitarity ai^e associated 
solely with the non-singular term Kj^\ and deviations from the universal radiation functions are as- 
sociated solely with the non-singulai" term A'^^^. In both cases, the generalization to higher orders is 
straightforward. 

In traditional parton showers, all the non-singular terms are dropped, and hence only the unitary 
singular stiucture remains. 



Exponentiating the leading singularities, we may replace them by the Sudakov factor. 



A(W,s,Qt) = exp 



(9) 



We thereby obtain the all-orders pure-shower Markov chain, 

Si{p}H,s,Ql,0) = A{{p}h,s,QI)5{0-0{{p}h)) 

H + exclusive above Q e 



fs d$t'l (10) 
E / . -ri^ Sr A{{p}H,s,Ql^,) Si{p}H+i,Ql+„Ql,0) . 
„ Jot d'P H 



H + I inclusive above Q e 



The shower may exponentiate the entire set of universal singular terms, or only a subset of them (for 
example, the terms leading in the number of colors Nc). More on the Markov formalism can be found 
in ref. [33]. We hope this brief introduction serves to put the developments below in context, and note 
that we will return to the restoration of the finite terms in the section on matching. 
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2.2 Dipole-Antenna Showers 



In leading-log dipole-antenna showers, the fundamental step is a Lorentz invariant 2 — 3 branching 
process by which two on-shell "parent" partons are replaced by three on-shell "daughter" partons. 
This 2 — )• 3 process makes use of three ingredients: 

1. An antenna function that captures the leading tree-level singularities of QCD matrix elements. 
This is the equivalent of the splitting functions used in traditional parton showers, with some 
important differences, as we discuss below. 

2. A kinematics map, specifying how the post-branching momenta are related to the pre-branching 
ones. This is the equivalent of the "recoil strategy" of traditional parton showers. 

3. An antenna phase space — an exact, momentum-conserving and Lorentz-invariant factorization 
of the pre- and post-branching phase spaces. Traditional parton showers, on the other hand, are 
based on a phase-space factorization which is only exact in the coUinear limit, and momentum 
conservation may only be imposed a posteriori. 

In the following paragraphs, we present the notation and normalization conventions that we shall use 
in the rest of the article for each of these pieces. 



Factorization: Labeling the three daughter partons i, j, and k, we write the integral over a three- 
body matrix element corresponding to that final state in factorized form as follows, 

m{p,,pj,p,)\ d$3 = \M2{S)\ d^2 [^^(^^p ' ^^^^ 

where the two-paiton matrix element we have introduced corresponds to the "pai^ent" configuration, 
in which we label the partons I and K. The branching process represented by this factorization is 
thus IK — )• ijk, with total Lorentz invariant Sijk = sjk = s. The IM3P/IM2P factor in eq. (11) 
represents the evolution kernel, whose (negative) exponential is the Sudakov form factor, cf. eq. (9). 

Phase Space: The dipole-antenna phase-space measure is thus [33] 

— — = dsij dsjk , = , (12) 

d$2 27r 167,2^^ {s,m],ml) 

where the Kallen function, 

A(s, mj,mj^) = + mj + mj^ — Ismj — 2sm\ — 2m'jm1^ , (13) 
in the denominator reduces to a/A = s for massless particles. 



Antenna Function: The ratio of matrix elements appearing in the integrand of eq. (11) is then 
decomposed into a symmetry factor, a coupling factor, a color factor, and an antenna function, 

\M3{pi,Pj,Pk)\^ 2 ^ -Q , s 

I^^^^^p = ^iK^ijk g tijk aij,,[s,Sij,Sjk) , (14) 



where S takes into account potential identical-paiticle factors as well as the possible presence of more 

ijk 



than one antenna in the parent (IK) configuration, is the relevant coupling factor, Cj,fc is a color 
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Figure 1: Contours of constant value of the antenna function, a^j/, for qq — )• qgq derived from Z decay 
as function of tiie two piiase-space invariants, with an arbitrary normalization and a logarithmic color 
scale. Larger values are shown in lighter shades. The (single) coUineai" divergences sit on the axes, 
while the (double) soft divergence sits at the origin. 



factor, and a^j^ is a generic color- and coupling-stripped dipole-antenna function, with superscript to 
denote a tree-level quantity. The thi^ee-particle matrix element is averaged azimuthally (over 0). Note 
that our use of lower-case letters for the antenna function is intended to signify that it corresponds to 
what is called a sub-antenna in ref. [36] for which lower-case letters are likewise used^. 

For illustration, contours of constant value of a^gg{s, Sqg, Sgg) as derived from Z decay are shown 
in fig. 1, over the 2 — 3 phase space, with an arbitrary normalization and a logarithmic color scale. 
This function is called in ref. [36] and is identical to the radiation function used for qq — qgq 
splittings in Ariadne. One clearly sees the lai^ge enhancements towards the edges of phase space, 
with a double pole (the overlap of two singularities, usually called soft and collinear) sitting at the 
origin, and single singularities (soft or collinear) localized on the axes. 

Writing the coupling factor as = Attos and combining it with the phase space factor, eq. (12), 
we have the following antenna function normalization 



That is, we use the notation a for the coupling- and color-stripped antenna function, and the notation 
a for the "dressed" antenna function, i.e., including its coupling, color, and phase-space prefactors. 

Note that (7^ x (phase-space normalization) leads to a factor as/ (47r) independently of the type of 
branching. As we believe that the formalism becomes more transparent if the origin of each factor 
is kept clear- throughout, we shall therefore use this factor for all branchings, instead of the more 
traditional convention of using Ug/ {2tt) for some branchings and a^/ (47r) for others. Obviously, this 
convention choice will be compensated by our conventions for the color factors and antenna-function 
normalizations, such that the final result remains independent of this choice. 

"Thus, in the notation of ref. [36], our dipole-antenna functions would be dj] — A''^, d}^ — d°, e° — 5-E3, /s ~ fi, and 




1 



(15) 



a 
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Color Factor: For color factors, we shall systematically adopt a convention optimized for shower 
applications, in which the color factor tends to Nc in the large- A'^c' limit whenever a new color line is 
created (rather than to Nc/2 as is the case for Cp in the standard normalization) and to unity when no 
new color lines are created (rather than 1/2 which is the case for Tr in the standard normalization). 
That is, we shall use 

Ca = Nc = 3, (16) 

Cf = ^r^ = 8/3 = 2CF, (17) 

Tr = 1 = 2Tr. (18) 

With all other normalizations fixed, the normalization of the antenna functions is now also unique. 
Indeed, with this choice, the radiation functions will turn out to be normalized in a way that makes 
their similarities more readily apparent. Another useful thing about this normalization is that the color 
factors now have a very simple interpretation. They provide a one-to-one count of the number of new 
color degrees of freedom that have been summed over in any given process. This gives a simpler 
counting and interpretation than in the standard normalization. Of course, the fact that there are only 
Nfj — 1 = 8 gluons, leads to corrections of order l/N^,, and here again it is trivial to let the "naive" 
color-line creation factor Ca be replaced by Cp for, e.g., qq — )• qgq, without artificially having to 
compensate by a factor of 2 in the radiation function — the eikonal part of the radiation function now 
remains invariant (as noted above in the discussion of the normalization of the radiation functions), 
and the difference in color factor is explicitly subleading, as it should be. A final argument is the 
question of which color factor to use, e.g., for a. qg ^ qgg emission. In the standard normalization, 
this could lead to confusion, since one parent would seem to "want" Cp and the other Ca, which 
differ by a more than a factor of 2. In the normalization used here, the difference between Cp and 
Ca is explicitly subleading in color, and hence it is clear that either of them could be used without 
any possibility for ambiguity at the leading-color (LC) level, again placing the proper difference in 
the proper place. 

To preempt confusion and illustrate how simple the translation between these convention choices 
is, consider the dipole-antenna function for gluon emission off a qq pair used in the Ariadne gener- 
ator [17]. As the symmetry factor is unity, this is just the matrix element squared for Z — 3 divided 
by the one for Z — )• 2 multiplied by the aforementioned phase-space factor, 

|M(Z^ggg-)P _ I2a, {l-y,,f + {l-yjk? ^ ^^^^ 



167r2s|M(Z s Svr yijUjk 
where yij = Sij / s = I — Xk- The factor 2as/ (Svr) in the first equation can be rewritten in two ways 

2as 4: as ^ ds 2as 8 as ^ cts 

— - = = Cp— or — - = = Cp— : (20) 

Svr 3 27r 27r Svr 3 47r An' 

it is purely for our own convenience that we choose the latter normalization. 

In a similar vein one could rewrite the DGLAP splitting kernels [37], which ai^e used in traditional 
parton showers [11, 12,38,39], as 

K-^rnj^^) ^ 1 as^^ l + z^ ^ las^^ 1 + ^21) 
Sij Sij2iT l — z s Air yjj(l — z)' 

^91-^9.9^^) ^ 1 as - z{l - z)r ^ la.^ ,,^ (l-^(l-^))^ ^ ^22) 

Sij Sij 27r z{l — z) s An yijz{l — z) ' 
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Global Dipole-Antenna (ARIADNE [17], GOG [36], WK [32], 
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Figure 2: Schematic overview of how the full coUinear singularity of parton / and the soft singularity 
of the IK pair, respectively, originate in different shower types. (Qj and 6^ represent angular vetos 
with respect to partons / and K, respectively, and Qjk represents a sector phase-space veto, see text.) 



where the gluon radiation function has absorbed a factor of 2 on the r.h.s. of the last line, due to the 
normalization choice. We note that, although these expressions look quite different from the dipole 
formula, eq. (19), they lead to identical singularities. This was shown in ref. [29] by identifying z as 
the Lorentz invariant energy fraction taken by the quark, z = Xi/{xi + Xk), and adding the radiation 
from the antiquark, qx QjQk- 

Shared Singularities: This examination of the different presentations of singularities brings us to 
the issue of "shared singularities". In traditional parton showers, as we have just seen, the full leading- 
log radiation pattern can only be obtained after summing over pairs of partons (which each radiate as 
independent monopoles), and care must be taken in the construction of the shower to make this sum 
approximately coherent to reproduce the correct singular behavior for soft wide-angle radiation. This 
dipole singularity is the simplest case of what we shall generally refer to as a shared — or multipole 
— singularity below; radiation whose full singularity structure (in a particular phase-space hmit) can 
only be recovered after summing over two or more radiators. 

A chain of such uniquely labeled and color ordered gluons, which could, e.g., represent a shower 
"event record" at a given point during its evolution, is illustrated in fig. 2. Below the schematic drawing 
we give an overview of how the full collinear singularity of parton /, and the full soft singularity of 
the IK pair, would be obtained for five different kinds of parton shower models, as follows. 

In a traditional parton shower, the full collinear singularity of each parton is contained in the 
DGLAP splitting kernel, P{z), that generates radiation off that parton. Since no other radiators share 
that collinear direction, there is no double counting at the LL level. (The kernel P{z) constitutes 
a complete subtraction term for the collinear singularities in real-emission contributions to an NLO 
calculation.) However, in this approach, the soft (eikonal) singularity between the IK pair must be 
obtained by summing the radiation functions of partons / and K together, and therefore it is essential 
in this type of approach that both the radiation functions and the shower phase-space factorization 
represent a correct partitioning of the soft region, with no so-called dead or double-counted zones. 

In the early eighties it was shown [40] that additional coherence effects can also be taken into 
account in this language, albeit approximately, by imposing angular ordering during shower evolu- 
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tion. This effectively represents a first step towards treating color-connected partons as multipoles 
in the shower context; partons I and K now effectively acquire some knowledge of each other, via 
their relative opening angle, and hence they no longer act as completely independent emitters. This 
improvement is denoted "Coherent Parton Showers" in the table in fig. 2. As indicated by the appear- 
ance of the function in the collinear term in such approaches, it is important to construct the angular 
ordering in such a way that the effect of the veto disappears in the collinear limit. 

In this paper, we follow the dipole-antenna approach to color coherence. This is motivated by 
the observation that, whereas the collinear singularities are associated with single logarithms, the 
parametrically larger double logarithms arise from soft (eikonal) factors. It therefore makes sense to 
change the underlying picture to a dual one in which parton pairs are the fundamental entities. Such 
pairs appear in the fixed-order literature under the name of antennce and in the shower one, under 
the name of dipoles. The latter term, however, usually means something else again in the fixed-order 
literature. To avoid confusion, we therefore call these pairs dipole-antennce . In this picture, the roles 
of soft and collinear singularities are interchanged, with respect to the parton picture. The soft double 
logarithms between neigboiing partons now come from a single term, which is thus guaranteed to be 
neither over- nor under-counted as no other pairs become doubly singular in the same phase-space 
region. The single logarithmic collinear radiation off a given gluon must now be partitioned among 
the two neighboring antennae that share it. (Note that quarks are still unambiguous in this picture.) 
The gluon case is represented by the line labeled "Global Dipole-Antenna" in the table in fig. 2. 

There is considerable freedom in how to partition the collinear radiation, because terms can be 
shuffled back and forth "across the gluon" while maintaining their sum constant. Two convenient 
examples are furnished by Ariadne [17] and by Gehrmann et al. (GGG) [36], which use different 
decompositions (see ref. [41] for some additional discussion of this point). 

The first important point concerns what to compai^e; obviously, the individual shower functions 
differ by collinear singular terms. Thus, if we naively subtract the Ariadne functions for qg — )• qgg 
and gg — ggg [17] from the corresponding GGG ones^ [36], we obtain 

^qg _ 1 ( 'iVik , yijVik , VikVik + y% ^ h ^ I \ l ( {l-y^jf + {l-yjk?\ 

"GGG-AR — 7 ^ — - \ — ^ - + -Vjk 



VijVjk Vjk Vij 2 2 J s \ VijVjk 




Vjk 

"^Vik , yijVik , yjkVik , 8^ 1 /^(i - yijf + (i - Vjkf 



(23) 



''ggg-ar — ~ ' I I I" 



yijyjk yjk yij s \^ yijyjk 



, yij + —-y^k + i^] , (24) 

s \ yij yjk 3y 

which differ by gluon-collineai^ singular terms. However, when we sum over the two possible order- 
ings of the gluons in eq. (23) and the three orderings in eq. (24), the discrepancies become coUinear- 
finite, 

^GGG-AR = "^GGG-AR + 0' ^ ^) = ~ ) (25) 
Ag'GG-AR = €'GG-AR+(i^^) + (^^j) = f • (26) 



^^Note that, in the shower context, different color orderings of the final state are represented as separate events, wherefore 
the shower function should only contain one of the permutations, corresponding to what GGG label a sub-antenna. 
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Sector populated by IK^ijk Sector populated by Jl^jki Sector populated by KJ->kij 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

yy = Sij/Sijk = 1 -Xfc yij = Sy/Sijk = 1 -X^ yij = Sy/Syk = 1 -Xfc 

Figure 3: The three phase-space sectors in a color-singlet gigjgk configuration, using p± as the dis- 
ciriminator for which sector a given emission/clustering/history belongs to. 



We see that, once summed over pemiutations, the Ariadne functions have substantially smaller 
finite terms than the GGG ones. The Ariadne shower is accordingly somewhat softer than the 
default ViNCiA one, which uses the GGG functions, but the singular terms are the same, in spite of 
the apparently singular differences between the individual shower functions. 

We shall usually choose the partitioning of GGG, which makes the collinear temis of the gluon 
antennae appear almost identical to those involving quarks in our parametrization, thus emphasizing 
their similarities. Note that for shower applications, the partitioning must be done in such a way that 
the resulting shower functions are positive definite. This is indeed the case for all the functions we 
consider in this paper — a counter example is given in ref. [43], where the positive-definite Ariadne 
antenna functions are, repartitioned a la GGG in a way that introduces negative regions in the individual 
99 ~^ 999 shower functions, while maintaining their sum constant. 

A different approach to the issue of how to partition the collinear singularities is to retain the 
full collinear singularity of the gluon in both of the neighbouring antennse, and combine this with 
phase-space vetos that allow only one or the other antenna to contribute to each given phase-space 
point, a possibility we have labeled "Sector Dipole-Antenna" in the table in fig. 2. The distinction 
between global and sector antennse is thus that in the former, several antennse (two that are singular, 
and possibly more that are non-singular) are summed over all of phase space, in a way such that 
their sum reproduces the full collinear singularity, and in the latter case every single term contains 
the full singularities, but only one term (the most singular) is allowed to contribute to any given 
phase-space point. An illustration of the sectors appropriate to one color ordering in the decay of a 
scalar, — )• ga9r9b^ is given in fig. 3. This approach has been suggested for use in NLO fixed- 
order calculations, and can be used with ViNCiA as well. Larkoski and Peskin (LP) [41] have also 
considered these kinds of antennae, including polarization effects. 

Finally, a different approach, which also treats dipole coherence exactly, consists of systematically 
partitioning both the soft and colUnear singularities of I and K into four terms, two of which treat 
parton / as the emitter, with K and H acting as spectator/recoiler, respectively, and the other two terms 
treating paiton K as the emitter, now with paitons / and L acting as spectator/recoiler, respectively. 
Catani and Seymour [25,44] labeled this a dipole model, but as this usage differs from an older use 
of the term dipole in parton-shower calculations to describe the sum of two such terms (in the context 
of the Lund dipole [15, 16]), we avoid confusion by referring to the CS type as a partitioned-dipole 
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shower and to the Lund-dipole/antenna type as a dipole-antenna shower. In the case of partitioned 
dipoles, the radiation off parton / is spht into two terms ("sides"), each of which contains "half" of 
the collineai- singularity of parton / and "half" of the soft singularity with either of the neighboring 
partons, e.g., 

ai^K = \co\\{I) + ^Soft(/K) , (27) 

where the subscripts are intended to denote that this is the term for I emitting and K recoiling. 
There are thus a total of 4 radiation functions involving partons / and K, but these terms can now 
all be constructed expUcitly so as to have the correct limiting behaviors. Obviously, there is some 
ambiguity concerning which functional form to choose for how to divide up the radiation among the 
various terms, which is why "half" is in quotation mai^ks. Recently, Nagy and Soper (NS) presented a 
proposal [45] for turning the CS subtraction scheme into a parton shower; several groups have since 
developed CS-style showers, most notably Dinsdale, Ternick, and Weinzierl [24] and the Sherpa 
group [23]. Although not based on the CS formalism, we note that the p^-ordered showers in Pythia 
8 are also closely related to this approach. 

Kinematics Map: The kinematics map specifies the details of how to reconstruct the parent mo- 
menta IK from the daughter momenta ijk and is equivalent to what is referred to as the the "recoil 
strategy" in parton shower language. In an old-fashioned parton shower [46], or a Catani-Seymour 
one [23-25,31,45], for instance, the recoil strategy usually implies classifying either I or K as the 
emitter and the other as the recoiler, with the recoiler being constrained to experience a momentum 
change only along its direction of motion in some frame (e.g., the rest frame of the emitter -i- recoiler), 
say, PiWpi- In dipole-antenna approaches, / and K can be allowed to share the emitter/recoiler roles 
more smoothly over the resolved parts of phase space, with a clear distinction only being made in the 
strictly singular limits. 

In principle, allowing recoils also outside the 2 — > 3 process itself, i.e., involving other partons 
than / and K, could be imagined, as long as the leading singular limits ai^e respected. This would 
change the subleading properties of the resulting shower approach, which might be deemed desirable 
in some contexts, although it of course would not alter the formal level of precision. However, it 
does make the fomialism more cumbersome, and hence we shall here restrict our attention to "local" 
recoil strategies, i.e., involving only the paitons / and K. With this restriction in mind, the constraints 
that must be fulfilled to obey the singular Umits and viable functional forms the kinematics map were 
discussed inrefs. [17,33,36,47,48] (including 2 — > n generalizations inrefs. [47,48]). 

For illustration, fig. 4 shows the branching phase space together with examples of the orientation 
of the post-branching partons in the CM of the branching antenna for various phase-space points, 
using an antenna-like kinematics map, the "Ariadne angle", according to which the two pai^ents 
share the transverse component of recoil. We shall return to kinematics maps later, but for the present 
merely note that the difference between an emitter-recoiler picture and the map used in fig. 4 is just 
an overall rotation (about an axis perpendicular to the paper), which vanishes in the soft and coUinear 
limits. 

3 A Shower Generator Based on Dipole-Antennae 

A parton shower algorithm can be constructed from knowledge of the Sudakov form factor A(Q^^, Q\2) 
representing the probability that no branching occurs between the scales Qei and Q%2- The Sudakov 
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Figure 4: Illustration of the branching phase space, eq. (12), for qq — )• qgq, with the original dipole- 
antenna oriented horizontally, an antenna-like kinematics map (the "Ariadne angle") in which the 
two parents share the transverse component of recoil, and (j) chosen such that the gluon is radiated 
upwards. 



form factor is in turn given by the exponential of a branching integral, 

A = exp(-^) , (28) 

where, following the conventions laid down in the previous section, the integral A of an antenna 
function over the 2 — )■ 3 antenna phase space is 

(We have suppressed the trivial integration over (/>.) Here, we have set all masses to zero, which 
approximation we adopt throughout this paper. PeifoiTning such integrals over all of phase space 
yields exactly the subtraction terms used in antenna-based fixed-order calculations, see, e.g., ref. [36]. 
(If the approach relies on phase-space vetos, such as in the case of sector antennse, we can treat these 
as step functions that are part of the antenna function a, so that eq. (29) remains valid). Within the 
shower algorithm, we need to evaluate the integral for a range of scales, and then invert the function 
to write the lower scale as a function of the Sudakov form factor. The two-dimensional nature of the 
integral means that we have to define coordinates, of which one will be the evolution scale of the 
shower. We can use this to our advantage, as we will see below, to allow for a number of different 
evolution variables within the same formalism. A functional inversion that is both flexible and efficient 
is accomplished by first using a simple overestimate of the antenna function, and then vetoing to obtain 
the exact result. 
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3.1 The Evolution Integral 



In eq. (29) we left the integration boundaries unspecified. For showering purposes, what is needed is 
not the integral over the entire phase space, but over a region bounded by two values of the evolution 
variable, 

A{s,Qei,Q%2)= I dsij dsjk a{s,Sij,Sjk) ; Qe2 < Qei ^ (30) 

which represents the integrated tree-level splitting probability when going ^om the scale Qei to the 
(lower) scale Qe2- This is the fundamental building block which we shall later exponentiate and 
invert in order to find the evolution equation, but first we need to recast it so that the evolution variable 
appears explicitly as an integration variable. 

We do this by a change of variables from the original invariants to a new set, one of which is 
the evolution variable, Qe, and the other we may label since it will play a role similar to, but not 
identical with, that of the z variable of traditional coUinear-based shower algorithms. That is, our 
generic evolution integral will have the form 

A{s,Q%i,Q%2)= I dQ E dC \ J\ a{s,Sij,Sjk) , (31) 

where \ J\ is the Jacobian associated with the transformation from {sij, sjk) to {Q'e, C)- 

One immediate difference between our approach and traditional collinear-based formalisms is that 
there is here no explicit dependence on the precise definition of it merely serves to (re)parametrize 
phase space, and from a set of generated (QIjiC)' the set of invariants {sij,Sjk) may be obtained 
without ambiguity. By contrast, in classic parton-shower approaches one would usually start from the 
coUinear-limit splitting functions P{z) or similar objects, and since these are only accompanied by 
an unambiguous definition of z in the collinear limit, the precise definition of z away from this limit 
(energy or light-cone momentum fraction? in which frame? with finite-momentum recoils?) results 
in an ambiguity which is not present in our treatment. 

There is of course a dependence on the functional form of Qe, which fomially enters starting 
from second order in the expansion for an IR safe obsei^vable (for IR sensitive ones, Qe needed 
as a regulator already at first order). Much effort has gone into debating and examining the vices 
and virtues of individual choices. Our stance is to order preferably in the inverse of the radiation 
function, since it is the singularities of this function which drive the logarithms; i.e., in p± for gluon 
emission and in the invariant uigq for gluon splitting to qq. Here, however, since we start directly from 
eq. (30), the phase-space factorization expressed in eq. (12) is preserved for any choice of Qe and 
and so rather than restrict ourselves to one specific form, we may instead seek a general solution to the 
evolution integral which will apply to a continuous class of Lorentz invariant evolution variables. By 
varying this unphysical pai^ameter (which effectively amounts to changing the factorization scheme 
since the evolution variable is what separates "resolved" from "unresolved" parts of the calculation at 
any given stage of the evolution) we obtain an estimate of the amount of scheme dependence which 
is generated by this variable — a dependence that higher-order matching will expUcitly reduce, as we 
shall see below. The remaining variation can then be interpreted as an uncertainty estimate. 

We stress that we here only give the possibility to vary this choice — further studies would be 
required to determine what a sensible range of variation would be, in the context of uncertainty esti- 
mates, in the same way that one discusses variations of other unphysical parameters for uncertainty 
estimates. The latter is obviously an art, not an exact science, but an art for which we can at least 
furnish the tools. 
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3.2 Trial Functions 



To simplify the equations, we shall make use of the veto algorithm to replace the integrand, a, by a 
simpler function, atrial, which we shall call the "trial function". Provided our trial function is larger 
than the actual integrand, the veto algorithm will allow us to recover the exact integral post facto, 
and if the overestimates are not extremely wild, only a small loss of efficiency should result. We will 
choose the trial function to have the same leading (double) singularities as the full antenna function. 
It will thus provide us with a simple Lorentz-invariant phase-space generator which is pre-weighted 
to take into account these leading double-logarithmic singularities of QCD; we may then reweight 
back down to more exact expressions efficiently. These more exact expressions may be complete 
matrix elements or dipole-antenna splitting functions, depending on what is available at a given order 
in perturbation theory. The overestimating function we shall use for all branchings is the following, 
with a normalization depending on whether the trial process is gluon emission or gluon splitting to qq, 



atrial-cmit = . "'^ Ca atrial-cmit = . Ca "( 99 QOO & C.C. (32) 



qgq 

Q99 
99 999 



as rf. - «s ^ Sijk I qg ^ qqq' & c.c. 

fltrial-split = ntlR atrial-spht = -. nflR < _ . (33) 

AiTSijk ^ 47rsijfc SijSjk [ 99 9m + QQ9 

where n j is the number of kinematically accessible flavours and Os is a parameter representing a 
"trial" a^; it should be greater than or equal to the latter. We will explain its use further in sec. 4. For 
comparison, the Eikonal (soft) approximation for gluon emission is 

q ( \ ^('^ ~ ^ij ~ ^jk) 'S ~ Sjj — Sjk _ , . ('XA\ 

i^eraitxS , Sij , SjP; ) — (^tiial-emityS, Sij , Sjj^ j , 

SijSjl^ S 

such that our approximation coincides with the Eikonal in the soft limit, Sij — and Sjk — )■ 0, 
and is lai^ger than the Eikonal everywhere else. Obviously, the normalization factor can be adjusted 
should extreme variations exceeding these upper bounds be deemed interesting to study; however 
we have not found that to be necessary in connection with any of the studies in this paper. Note 
also that using the Eikonal for gluon splitting is likely to give a very large overestimate over most of 
phase space, as compared to the physical antenna functions (which only generate single logs, while 
the Eikonal generates double logs), thus leading to low efficiency in the subsequent veto step. This 
may therefore be a technical point to improve on in future work, but for the time being we prefer the 
simplicity of having the same functional form to work with for all trials. Let us emphasize that the 
veto algorithm ensures that there is no trace of the overestimator present in the final results, either in 
tree-level expansions or in the Sudakov exponentials. The only sensitivity to the overestimator is in 
the speed of the calculation, which we have tested to be comparable to that of standalone Pythia 8. 

Inserting the trial function atriai-emit from eq. (32) and the Jacobian for the transformation from 
{sij,Sjk) to arbitrary (Ql-, C)> eq- (30) then becomes 



Arial-cmit(s,Q|l,Q|2) = / "^dQldClJl ^""^ , (35) 

4vr SijSjk 



where we have kept Og inside the integral since the renormalization scale may vary over phase space. 
The last remaining step before we can solve this equation is now to rewrite the term \ J\/{sijSjk) in 
terms of Qe and ( for a sufficiently general class of functional forms of Qe- 
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C definition C definition (Log) 




0.0 0.2 0.4 0.6 0.8 1.0 -4 _3 _2 -1 1 2 

yy = Sij/SijK = 1 -Xi, Ln(yi2) = Ln(Si2/Si23) = Ln(1 -X3) 



Figure 5: Illustration of the ( definition, eq. (36). The physical phase space, shown in grey, is the 
same on both panes, and the ( definition is also the same, but on the left the phase space is shown on 
a linear scale in the branching invariants and on the right on a logarithmic one. 



We shall start by specifying (. In principle, one could define a different ( for each possible choice 
of Qe, but this would be cumbersome and would not lead to any significant gains. We have therefore 
settled on one particular form for ( to use for all evolution variables in this paper. There aie two 
requirements for such a ( definition: 

• It should be linearly independent of any Qe that we would conceivably consider. I.e., C should 
not itself be a candidate for an evolution variable. 

• Curves of constant ( should intersect curves of constant Q e once and only once for all > 0, 
such that the Jacobian is well-defined. 



A C definition that fulfills these conditions for the entire class of Qe variables we shall consider 
is the following simple ratio of invariants, illustrated in fig. 5: 

Q = ^^ ^ l_C = ^j^_. (36) 

Sjj + Sj}^ Sij + Sjfc 

We emphasize again that all we have done so far is recast the Lorentz invariant phase space in terms of 
two new variables which are, themselves, arbitrary. There is no explicit dependence on the particular- 
form of (. (There is a dependence on Qe; of course, but only through the boundaries of the integral.) 
To compute the Jacobian, we will need the derivatives, 

dC_ ^ s,k ^ C(l-C) _dC_ ^ -s^j ^ -C(l-C) ^37^ 



The Jacobian is. 



det 



d{sij , Sjfc} 

d{QhC} 



det 



\d{s,j,s,k}) ({1-0 



dQ' 



E , „ dQl 



ds 



jk 



(38) 
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Inserting this in eq. (35) we get a master equation for evolution in an arbitrary variable Q e (for 
trials distributed according to the function atriai-cmit): 



Ca 
47r 



2 r, 



Cn 



2a, 



ax(Q|;) 

.{Ql) C(l 



C) 



Si 



dQl , dQl 



n -1 



+ Sjk 



ds-i 



jk 



(39) 



where the function in square brackets represents the leftovers from the Jacobian, and the functions 
Cmin(Q|;) ^nd Cmax(Q|;) re-express the phase-space triangle in terms of Qe and (. Given any specific 
form of Q^; these three functions can be derived, and hence a particular evolution equation obtained. 



3.3 Evolution Variables 

Given the structure of eq. (39), one sees that the evolution integral will become particularly simple for 
any evolution variable which satisfies the following simple differential equation. 



Qlis, Sij, Sjk) — Ke 



\..dQl 



+ s 



jk 



dsjk 



(40) 



Rather than base our formalism on one particular choice of evolution variable, as is usually done, we 
shall therefore instead derive our formalism so that it applies to any evolution variable which satisfies 
eq. (40). Making only this requirement, the evolution integral simplifies to 



-^trial-cmitls, QIi: QI2) 



KE 



Ca 



?2 



2 
E 



Q 



Cmax{Q|;) 



Cn 



dC 



2a, 



C) ' 



(41) 



which may be solved once and for all, independently of the choice of Q^;. (The explicit dependence 
on the form oIQe will reemerge, as it should, when translating from a generated set of {Qe, C) back 
to the branching invariants, (sjj, Sjk), but this is a separate step, to be treated below.) As explained 
later, we will take a, to depend solely on Ql, and to be independent of C,. 

Though any symmetric power series in the branching invariants is a solution to eq. (40), for the 
purpose of this paper we shall impose two additional "reasonable" boundary conditions. Firstly, that 
an infinitely soft branching should always be classified as "unresolved" for any finite value of the 
evolution variable, i.e., the evolution variable must go to zero when both invariants vanish, Qe{^, 0) = 
0. We note that this apparently mild restriction will nonetheless prevent us from considering a specific 
class of variables, which in fact include angular ordering, since in an angular-ordered cascade, a large- 
angle emission may be resolved, even when soft. We shall also require that the evolution variable be 
symmetric in the two invariants, QE{yij,yjk) = QE{yjk,yij)- 

We can put the differential equation (40) in dimensionless form by dividing by s. 



yE{yij,yjk) = ke\ — h 



where Ye = Ql/s, iHj = Sij/s, and yjk = sjk/s. 
If we define 

L+ = (In yi + In y2)/2 , L_ = (In yi - In y2)/2 , 
then eq. (42) takes the form. 



Ye(L+,L„) = KE 



dVE 



(42) 

(43) 
(44) 
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ViNCIA 



a=1/4,p=1 



a=1/2, p=1 



a=1, p=1 



a=2,p=1/2 




0.0 0.2 0.4 0.6 0. 

yii = Sj/Sin = 1 -X, 
a) V-ordering 



yii = Sj/Sin = 1 -X, 
b) mo -ordering 



0.0 0.2 0.4 0.6 0.8 1.0 



c)px -ordering 



0.0 0.2 0.4 0.6 0.8 1.0 

yii = Sij/Siik = 1 
d) iJJ-ordering 



Figure 6: Illustration of the a parameter in eq. (47). Note that the second plot from the left corresponds 
to choosing the daughter antenna mass as the evolution variable, and the third pane corresponds to the 
Ariadne definition ofp±. Contour labels indicate values of ye = Q\/s. 



which has the general solution, 

YE{Sij,Sjk) 



/o(L_)exp(L+/K£;) + Co 

/o(s.i/s,fc)(Si,S,fc/s2)V(2Ks)+Co. 



(45) 



where /o is a dimensionless function satisfying /o(x) = /o(l/x) and the requirement that Y£;(0, 0) = 
0, and Co is a constant (which we uniformly set to zero). 

For illustrative purposes, it will be convenient to introduce additional parameters — a main shape 
parameter, a, and two auxiliary parameters h and p — and take 




a\ p 



(46) 



This yields a continuously deformable class of evolution variables that fulfill eq. (40) subject to our 
two additional conditions; the corresponding expression for the evolution variable is. 



+ Sjk 



Sjk 



a\ p 



subject to the constraints. 



a> , b <l , p> , 2ap 



KE 



(47) 



(48) 



The overall normalization is fixed so that the maximum of the evolution variable is the invariant mass 
of the dipole: a sets the relative soft/collinear resolution power (to be explained further below); b ^ 1 
allows for variables which do not go to zero on the axes (i.e., for which purely collinear branchings 
may appear resolved, as is the case, e.g., for energy ordering), and thus requires an additional infrared 
regulator independent of the evolution variable; and p allows modifying the overall speed of the 
evolution over phase space. 

The effect of varying a is illustrated in fig. 6, with 6 = 1. We increase a from left to right in the 
figure. Small values, toward the left, yield evolution variables that are "better" at resolving phase- 
space points towards the origin (corresponding to soft branchings) whereas large values of a yield 
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ViNCIA 



a=1, b=0, p=0.5 a=1,b=0, p=1 a=1 , b=0, p=2.0 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

yai = Sar/Sa,b = 1 -Jffi Xa, = SjS^^b = 1 -Jffi Ya, = Sa,/Sarb = 1 

1 ^ 2 

e) E '-ordering e) iJ' -ordering e) E -ordering 



Figure 7: Illustration of the b and p pai^ameters in eq. (47). Note that all of these variables are propor- 
tional to the energy of the emitted gluon (in the CM of the dipole-antenna) to some power. Contours 
indicate values ofi/E = Q'e/s. 



variables that are "better" at resolving points neai^ the axes (corresponding to purely coUineai^ ones). 
A shower based on the evolution variable to the far left in the figure would generate soft branchings 
earlier in the shower than coUinear ones, while a complementary ordering would result in a shower 
based on the variable illustrated on the fai^ right. We stress that the leading limiting behavior in the soft 
and collinear regions are in all cases the same, but these differences will lead to subleading differences 
between the showers. The size of these differences may be estimated by varying the evolution variable 
used to generate the showers. 

An illusti'ation of the b parameter is obtained by comparing the plots in fig. 6 to those in fig. 7, 
where the former all have 6 = 1 and the latter all 6 = 0. As mentioned above, choosing b ^ 1 gives 
evolution variables that do not go to zero on the axes. Accordingly, the contours on the plots in fig. 7 
intersect the axes, while the ones in fig. 6 do not. The corresponding showers would generally have 
greater sensitivity to the infrared region and thus to hadronization effects. These variables correspond 
to those that are used in some analytic resummation calculations [49]. Finally, going from left to 
right in fig. 7, we see that small values of the speed parameter, p, correspond to a faster progress of the 
evolution variable over phase space, whereas large p values give the opposite. This speed has no effect 
on the generation of the first branching, but it does affect the value of the restart scale for subsequent 
branchings, which we will return to below. 

Some specific examples of evolution variables of this form that could be useful for Monte Carlo 
purposes are given in table 1, where we also give the corresponding C, limit functions and show how 
the generic functional fonii, eq. (47), simplifies considerably in several cases. For instance, setting 
a = l, 6=1, p = l gives the p^-ordering variable used both in Ariadne as well as in more 
recent work [32,33]. Contours of constant value of this variable are shown on the third pane of fig. 6. 
Contours illustrating the other variables in tab. 1 can also be found in figs. 6 and 7. (£'|,^-ordering is 
shown for n = 1 only, for which case we leave out the explicit subscript n.) Note that the name E^^ 
does not imply that this variable represents a physical transverse energy, but rather that it represents 
an interpolation between p_\_ and E*, with lower values of the n parameter making it closer to and 
higher values making it closer to E* . The V variable is named simply for the shape it has over phase 
space, like a V pointing towards the soft region, cf. the leftmost pane in fig. 6. 

For completeness, a few important examples of evolution variables that are not covered using our 
formalism are the traditional 1 — )• 2 parton-shower ones in Pythia and Herwig, illustrated in fig. 8. 
For the forms of the HERWIG and Pythia evolution variables, translated to our phase-space notation. 
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Name 


a 


b 




Resulting form for Q'^ 


Siniii \ ^ E J Sina-x E J 


1 


-ordering 


1 


1 


1 


s 


'^^'XAriadno 


lT^/l~Ql/s 
2 


2 


m/) -ordering 


1 

2 


1 


1 


2mm{sij,Sjk) 


= 2m% 


^ E 1 
2s ^ 2s 


3 


i?* -ordering 


1 





1 


S 


= AEf 


1 


4 


F-ordering 


1 
4 


1 


1 


\J s{Sij + Sjk) — 


\J s\Sij — Sjk\ 


2 


5 


£'|,„-ordering (n > 1) 


2n 


1 


1 

2n 


n = I : 


s 


2 



Table 1 : Examples of evolution variables in the fonii of eq. (47) and corresponding to the illustrations 
in figs. 6 and 7. The nominal C, boundaries for E* ordering would lead to infinities, so for practical 
applications the bounds implied by the hadronization cutoff should be used instead. 



Jetset & Fortran Pythia Pythia 6.3-i- & Pythia 8 Herwig-i-i- 



Virtuality-Ordering: side a p^e„o|-Ordering: side a Angular Ordering 




yar - Sar/Sarb - ^ -^b /ar - Sar/Sart} - "1 -^b /ar - Sar/Sarb - 1 -><b 



Figure 8: In an old-fashioned parton shower or partioned-dipole shower, the process ah — )• arh is 
divided onto two terms, one representing emission off parton a and the other emission off parton b. 
Left Pane: contours of constant Sar, i-c, the virtuality that corresponds to a* ar in a virtuality- 
ordered parton shower. The inset shows the equivalent contours for emission off side h. Middle 
Pane: contours of constant pxcmh the variable used in the p^-ordered Pythia shower. Note that 
for the virtuality-ordered shower, additional vetos on the emission angle, not shown here, must be 
imposed to enforce coherence, while in the -ordered case, this is less crucial due to the use of dipole 
kinematics. In an angular-orA&r&A parton shower (right pane), each parton is still evolved separately, 
but the potential for double counting has been removed by effectively restricting the emission from 
each parton to non-overlapping regions, here angular-ordered cones, and hence we can represent the 
two terms on one and the same plot. (Note: while the original Herwig implementation of angular 
ordering did imply some overlap in the soft region, this has been removed in Herwig-i-i-.) The price 
to pay is that this introduces an artificially unpopulated dead zone in the phase space, illustrated by 
the striped area. The contour labels denote values of ue = Q%/s. 
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Prsym-Ordering m^msym = d -T)-Orclering 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

/ar - Sar/Sarb - 1 -^b /ar - Sar/Sarb - 1 -^b 



Figure 9: Examples of evolution variables for which our C, definition would be multi- valued for each 
value of Qe- These particular ones were obtained by merely symmetrizing p_\_ (left) and m/j (right) 
in all three branching invariants. 



we used the following (for evolution with parton I as the emitter): 

Pythia : m|, = Sij 

virtuality-ordering 

PYTHIA : 4pL,„i, = 4^(«-«ifc)(^ii + «J-fc) 



{s + SijY (49) 



Pxevoi-ordering 

HERWIG++ : ql i = 4s ( — ■ r ) = 4s 



2 / X 2 



angular-ordering 



where we used Xi = {1 — sjk/s) in the last expression. For the Pythia p_lcvo1 variable, j is the 
emitted parton, and k is the recoiling one. Analogous expressions hold when parton K is the emitter, 
with the substitutions Sij -f-)- sjk- 

Likewise, variables that do not have well-defined Jacobians with our C, definition, such as those in 
fig. 9, cannot be used without additional work. Examples include, 

Pisym = 27^^^ (50) 



m 



2 

DSym 



3mm{sij,Sjk,Sik) = 3s(l - T3) , (51) 



where T3 is the Thrust event shape variable for a three-parton configuration. The p±sym variable is 
similar to the C-parameter for a 3-parton configuration, 

Cs = 6- JllIlhfUL. ^ < 3 (52) 

[s - Sij){s - Sjk){s - Ski) 4 

which one could therefore also imagine using as the basis for an evolution variable, with = IsCs. 



3.4 The C integral 

In general, one might expect the running coupling to depend on ( as well as Q^; and s. If we 
simpUfy the dependence by taking it to depend only on Q^; and/or s, then the integral over ( is 
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ViNCiA Evolution Windows 
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Table 2: The evolution windows used in ViNCiA, with the Qe boundaries and active number of 
flavors corresponding to each. The number of active flavors is the same for windows 2 and 3, but the 
boundaries for trials are different, due to the different Q Emm values. This improves the efficiency 
of the generator. The first window will not actually extend down to zero in practice, but will instead 
be cut off by the hadronization scale. 



straightforward: 

/•Cmax(Q|) 1 

^c(Cmin(Q|), Cmax(gl)) = / dC — -- = In 

One would expect any ( dependence in as to be reduced by computing to higher orders in perturbation 
theory. At a fixed order, a (^-dependent as could still be accommodated, for example, by choosing a 
(^-independent (large) for trials and then applying the veto algorithm. The case of a Q £;-dependent 
Us will be treated explicitly below. 

3.5 Evolution Windows 

To simplify the trial generation further we shall generate trial branchings in a larger phase-space 
region than the physically allowed one, again using the veto algorithm to avoid generating any actual 
branchings in the unphysical region. Specifically, we divide the generation of trial branchings into 
discrete windows in Qe — given in table 2 — and, in each such window, replace the C limits in the 
previous equation by constant ones, 

CminiQE) ~ Cmin(Q£;min) ' Cma,x{Q e) ~ CmaxCQiJmin) i^'^) 

where QEmin is the value of Qe; at the end of the current window (e.g., the next flavor threshold or, 
ultimately, the hadronization scale). If none of the generated trials fall within the current evolution 
window, the evolution should be restarted at Q^; = QEmin, upon which the QEmin and ( boundaries 
should be updated to correspond to those of the next evolution window. If the crossing corresponds 
to a flavor threshold, the normalization of the trial function for gluon splitting should be updated with 
the new number of active flavors, and if depends explicitly on Qe, then Aqcd should, likewise, be 
updated. 

3.6 Qe integral for QE-independent 

With the C integral for trial branchings having now effectively become a constant depending only on 
the current "evolution window" (i.e., the current Qmin). we may perform the Qe integration indepen- 
dently of the ( one. We do this first for the case where the renormalization scale in is constant over 
the branching phase space. 



Cmax(Q|;)(l Cmin(Q|;)) 
Cmin(Q|;)(l ~ Cmax(Q|;)) 



(53) 
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In this case, the evolution integral, eq. (41), becomes 



-4trial-cmit(Qsi) Q£;2) = 2asK£;^^/^(Cminj Cmax) In I I ' (^5) 

where Cmin.max are the ones appropriate to the current evolution window, as given in tabs. 1 & 2. 

The equivalent expression for the trial function for gluon splitting, ^triai-spiit> only differs by an 
overall normalization factor, cf. eqs. (32) & (33), and hence we do not reproduce it explicitly here. 

In addition to a zeroth order (fixed) as or a running q^j that depends only on the parent dipole- 
antenna mass, as{s), these expression will be needed, for instance, to use as{p±) together with any 
non-p± evolution variable. Technically, we accomplish this by setting the trial Og equal to unity (or 
some other relatively large number) in the above equation and then accepting the generated {QeX) 
pair with the probability as{pA_)/dis- 

3.7 Qe integral for first-order QE-dependent 

We shall also allow for the possibility to use a first-order running as, with a renormahzation scale that 
depends explicitly owQe, 

"^^^'^^^ = 6oln(A;2Q|/A2) ' ^^^^ 

where 

^ nC. - 2n,fn 

" 127r ' 

and where is an overall scaling factor, n j is the active number of flavors, and A is the appropriate 
(n/ -dependent) value of Aqcd- 

The evolution integral then becomes 

A (,n^ ^ - C-A 1 fQli dQl 1 

Ca 1 

= 2ke—. — ItyCrmn-iQm&xjT- / 



<lx\ 1 
x| ln(x|) 

"si) 




where 



2K^^/c(Cmin, Cmax)^ In f ) , (58) 



4 = ^ . (59) 



3.8 Qe integral for QE-dependent 

More generally, we can incorporate a Q^; -dependent as by changing variables using the /3 function, 

dasiQl:) 



dlnQl 
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/3(a,) . (60) 



This allows us to rewrite eq. (41) as follows, 



2 „2 N Ca /■°=(Q1;2) das /"CmaxCQlK)) 



When we make the substitutions of eq. (54) here, the inner integral will be independent of a^. With 
the two-loop beta function, 

P{as) = -boal - bial , (62) 

where 60 is given by eq. (57) and 

= 24^^^ ' ^^^^ 

the Os integral is simple, 

MQl,) das _ l^J asiQh) \ I I ^J ^ + bi/boasiQh) \ (^4^ 



"s{Q|i) bo \as{Q%2)' ^0 \l + bi/boas{Q%i) 

This function can be inverted readily using a Newton-Raphson solver, which can likewise be used to 
obtain Q\{as)- It can be extended readily to higher loops because additional orders only introduce 
new denominator factors of the form (1 + ca^). 



3.9 The Evolution Equation 

We now have all the pieces in hand to construct the evolution equation for a generic shower sub- 
ject only to the conditions outlined in the preceding paragraphs; that the individual trial branchings 
be 2 — )• 3 mappings from on-shell momenta to on-shell momenta, respecting the Lorentz-invariant 
phase-space decomposition, eq. (12), for any evolution variable that satisfies the differential equation, 
eq. (40). 

The generating function for such a shower is the Sudakov form factor, 

A(g|i, gL) = exp (-^(Qli, QI2)) , (65) 

where we may substitute for A either of the expressions eqs. (55) or (58). 

In order to generate trial branchings according to this Sudakov, we must solve the equation 

i? = A(Q|i,Q|2) (66) 

for Qe2, where is a random number distributed unifomily between zero and one and Qei is 
the "(re)starting scale". The latter represents the scale the shower is being restarted at prior to the 
generation of the next trial branching. To give an idea, this can either be the full dipole center-of-mass 
energy, y/s, which will usually be the case for the very first branching following a resonance decay, 
or, later on in the shower evolution, the scale of the preceding trial branching. 

The solution of this equation is the paramount reason we chose to use a simplified antenna func- 
tion for trial generations. It would have been possible to solve the evolution integral itself for more 
complicated trial functions, but the inversion of eq. (66) to solve for Qe2^^^ function of R and Qei 
would then have been much more cumbersome. 

Solving the evolution equation for a -independent a^, using eq. (55), yields 

QI2 = QliR' , (67) 
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with 



47r 



2dsAvsCA/c(Cmin(Q|;mm)' Cmax(<3|;mm)) 



(68) 



Solving the evolution for a first-order running Q^; -dependent as, using eq. (58), yields 



with 



b' 



^2 _ ( ^JiQeI 



UCa - 2nfTR 



2K£;CA/c(Cmm(Qsmin)' Cmax(Q£;min)) 127r 



(69) 



(70) 



As mentioned earlier, given any set of branching variables {Q%, C) we may obtain the invariants 
{sij, Sjk) without ambiguity. Thus, the next step is to generate a ( value, given Qe2- Since we defined 
C independently of Qe, we may do this once and for all, with the solution applicable to all evolution 
variables. To generate a random ^ value distributed according to the integrand of the integral, 
eq. (53), we must solve the equation 



Rt 



^( ( Cmin ) Cmax ) 



(71) 



for where i?^ is another random number uniformly distributed between zero and one, the integral 
given by eq. (53), and Cmin(Q£;min) is given by the evolution windows (Table 2) and by the evolution- 
variable-dependent ( limits (Table 1). 

We solve eq. (71) by first translating to the variable r. 



Cn 



Cn 



J- ^max Sjmin 



generating a random value for r 



and finally solving for C, 



r = r„ 



r \ 



1 + r 



(72) 

(73) 
(74) 



If the ( generated in this way falls outside the physical phase space, 

C<Cmm(Q|) V C>Cmax(Ql;) (75) 



the branching is vetoed. This occurs some fraction of the time for the simple reason that we generated 
the trial branchings in a hull encompassing the physical phase space. That is, the trials ai^e generated 
on a phase-space region bounded by Cmin,max(Q|;min)' whereas the physical phase space at Qe is 
bounded by Cmin,max(<5|;)- Since the physical branching probability outside the physical phase space 
is obviously zero, the probability to accept unphysical trial branchings should be zero as well. This is 
accomplished by the veto. 

The last step is to obtain values for the pair of phase-space invariants [sij, sjk) in terms of which 
we cast the original evolution equation, eq. (30). Since different forms of Qe depend in a different 
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way on these invariants, this step is obviously (5£;-dependent. Here, we give the inversions relevant 
to the four evolution variables so far implemented in ViNClA. 



with 



Sij — Cffevolution(s, C) j Sjfc — (1 ~ C)5evolution('S, Q^;, C) ("76) 



,,, ^ Qe\/s 

Typel: g„, = — , ^ . (77) 

Type 2: g^^ = —-9^-— (78) 

mm(C, 1 - C) 

Type 3: qe* = QeV^ , (79) 

Type 4: gy = — , . (80) 

2.(min(C,l-C) + ^T^) 

More generally, in terms of the function /o that parametrizes the general solution (45), the inversion 
takes the following fonii. 



s I 



2 \ KB 



3.10 The LL Shower 

In the previous subsections, the ingredients for generating a single trial branching with a trial branch- 
ing function atrial were described. To obtain an LL shower, it suffices to accept each trial branching 
with a probability 

— —-^ — r , (82) 

O^s Lijk Q^triaU'S, Sjj, Sjk) 

where the a^/as ratio takes into account the possibility that the trial generator could be using a 
nominally larger than the physically desired one, the C/C factor represents the same for color 
factors, and the antenna function ratio matches the trial function onto the desired physical splitting 
antenna for the relevant 2 — >• 3 branching. Because we chose atrial to be larger than (or equal to) the 
true antenna function Sll everywhere in the dipole-antenna phase space, this probability is always 
less than (or equal to) 1. We must also require oll to be non-negative in order that the ratio here be 
interpretable as probability. 

We initiate the shower in electron-positron collisions with quark-pair creation from the interme- 
diate vector boson. At each stage in the shower, a gluon will be emitted, or a gluon will split into 
a quark-antiquark pair. The shower itself evolves in the leading-color approximation, so after each 
emission, the number of different antennse grows by one, whereas each splitting leaves the number 
of antennse unchanged. We must allow all the different antennse to branch, of course; we do this by 
computing trial branchings for all of them, and picking the antenna with the highest trial branching 
scale. For that antenna, we then apply a veto with the probability given in eq. (82). 

When a branching is accepted, the physical replacement of partons I and K by i, j, and k in the 
event record next has to be performed. It is here that the dependence on the kinematics map enters. 
Our treatment of this point is identical to that described in detail in ref. [33], and the implementation 
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in ViNCiA retains the possibility to choose between the three different maps defined there. These 
maps all have identical limits in the LL singular regions, but differ from each other elsewhere. 

For a pure strongly-ordered LL shower (i.e., a level of approximation comparable to all other 
currently existing shower Monte Carlo implementations), the evolution should then be restarted from 
the scale of the current trial (regardless of whether that trial was accepted or not, as per the veto 
algorithm; the phase space above that scale has already been probed, and hence — according to the 
strong-ordering condition — should not be probed again.) 

To examine the quality of this type of approximation independently of the shower generator, we 
use Rambo [50] (an implementation of which has been included in ViNCiA) to generate a large 
number of evenly distributed 4-, 5-, and 6-parton phase-space points. For each phase-space point, 
we evaluate the Z — n matrix element using MadGraph (suitably modified to be able to switch 
subleading color terms on and off). We then compute the tree-level LL antenna-shower approximation 
corresponding to the same phase-space point, based on nested products of 2 — 3 branchings subjected 
to the condition of ordering in the chosen evolution variable. Finally, we form the ratio between this 
approximation and the full matrix element. This is similar to what was done in ref. [34]; where 
that study was limited to the emission of two partons, the addition of a new automated interface to 
MadGraph allows us here to extend the corresponding comparisons through four emissions, thus 
making it possible to illustrate in detail how the agreement or disagreement changes with increasing 
number of emissions. 

Using a clustering algorithm that contains the exact inverses of the ViNCiA 2 — 3 kinematics 
maps [34], we may perform m clusterings of the type k) — (/, K) in a way that exactly re- 
constructs the intermediate (n — m)-parton configurations that would have been part of the shower 
history for each ?i-parton test configuration, for any of the three kinematics maps so far implemented 
in ViNClA. This gives us an exact reconstruction of how the antenna shower would have populated 
each path. The strong ordering condition corresponds to step functions in the shower approximation. 
E.g., for Z qig29mi, we have [34], 



R. 



e(Q3A - QiA)aqg{l, 2, 2,)agg{l2, 23, 4) + @{Q^b - QiB)agq{2, 3, 4)a<j,-(l, 23, 34) j \M2{s) 

|M4(1,2,3,4)|2 

(83) 

where hatted variables fj denote clustered momenta, a denote 2 — 3 antenna functions, | M„ | ^ denote 
the color-ordered ?7,-parton matrix elements, s is the total invariant mass squai^ed of the n-paiton 
system, and the ordering conditions depend on 

Qaa = Qe(1,2,3) ; Q3A = Qi<;(12,23,4) 
QiB = Qi?(2,3,4) ; Q3B = Qi?(l,23,34) 

The numerator of eq. (83) thus reproduces the shower approximation expanded to tree level, phase- 
space point by phase-space point, for an arbitrary choice of kinematics map, (i, j, k) — (rj,jA;), and 
evolution variable, Qe- R thus gives a direct tree-level measure of the amount of over- or under- 
counting by the shower approximation, with values greater than unity corresponding to over-counting 
and vice versa. 

We use the kinematics maps defined in ref. [33], 

^4 + ^l. 




v^ps = <! _ \ :!'^:<'^:'' (86) 
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Figure 10: Strongly ordered parton showers compared to matrix elements. Distribution of 
logj^Q (PS/ME) in a flat phase-space scan. Contents normalized by the number of generated points. 
Spikes on the far left represent the underflow bin — dead zones in the shower approximations. Gluon 
emission only. Matrix-element weights from MadGraph [51,52], leading color (no sum over color 
permutations). 



where On^ is the angle between the after-branching parents in the CM frame of the branching. We 
show the results of these comparisons in fig. 10, for four different shower approximations: 

• GGG: -ordering using default ViNCiA settings, i.e., the GGG antenna functions and the 
ipAR kinematics map for all branchings. I.e., the parents share the recoil in proportion to their 
energies in the CM of the dipole-antenna. 

• V'PS p_L-ordering using the GGG antenna functions and the parton-shower-like (PS) longitudinal 
kinematics map. I.e., the parent with the largest invariant mass with respect to the emitted parton 
recoils only longitudinally. 

• m£,-ord: m/j-ordering using the GGG antenna functions and the tpAR kinematics map. 

• ARI: -ordering using our best imitation of the what the real Ariadne program does. It uses 

-ordering, but with the Ariadne radiation functions instead of the GGG ones, and it also 
uses a special recoil strategy, as follows; for qg dipoles, the quark always takes the entire recoil 
(in the CM of the dipole), whereas for gg dipoles, the V'ar angle is used to distribute the recoil. 

In aU cases, we compare to one leading-color (LC) matrix element, i.e., before summing over colors, 
and with all color factors having been divided out. We present an extensive set of comparisons for 
different ordering vaiiables in appendix A. 

The two bins around zero correspond to the parton-shower approximation having less than a 10% 
deviation from the full matrix element. At four partons, on the left-hand pane, these two bins contain 
over 35-60% of the sampled phase-space points, depending on the approximation, with tails extending 
out towards larger deviations. The spikes at the extreme left edge of the plots represent the underflow 
bin, including — oo, which coixesponds to zones in which all of the possible shower histories have 
been removed by the strong-ordering condition. Such dead zones are characteristic of (ordered) LL 
parton showers, when the ordering variable is more restrictive than pure phase space. We shall later 
discuss how to remove them while simultaneously improving the approximation in the ordered region 
as well. 

For all multiplicities, the default -ordering with the antenna-like Ariadne recoil map appears 
to generate the best overall agreement (narrowest distribution). The parton-shower-like longitudinal 
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Figure 1 1 : Strongly ordered parton showers, using the A RIADNE radiation functions and two different 
recoil strategies, compai^ed to matrix elements. Distribution of log^g (PS/ME) in a flat phase-space 
scan. Contents normalized by the number of generated points. Gluon emission only. Matrix-element 
weights from MadGraph [51,52], leading color (no sum over color permutations). 



recoil map (thin solid line labeled V'ps). following the spirit of Pythia 6 and showers based on CS 
partitioned dipoles and the dipole-mass ordering (dashed line labeled mi^-ord) give slightly worse 
agreement (wider distributions). 

The "ARI" case (thick solid line) has no dead zone for this process (due to the special kinematics 
map), but it also appears to generate a somewhat wider, and systematically softer (shifted to the left) 
distribution, than the GOG ones. To examine further whether this is an effect of the intrinsically softer 
radiation functions used in Ariadne (as shown in ref. [43]), or of the special recoil strategy employed 
by it, we plot in fig. 1 1 the result for 4 and 5 partons using the Ariadne radiation functions, either 
with the default Ariadne recoil strategy (using ipAR for gg dipoles only, with the quark recoiling 
for qg ones — thick solid line, as above) or just using ipAR for all dipoles (solid histogram). From 
this comparison, we conclude that there is no significant advantage in assigning all the recoil to the 
quark in qg dipoles. Although it removes the dead zone, together with a tail of largely undercounted 
events, the peak is actually degraded slightly, which we interpret as indicating a worse agreement with 
the matrix elements close to the singular regions. Physically, this could be consistent with the quark 
direction having to remain unchanged when the gluon branches coUinearly. The position of the peak, 
slightly to the left of zero, however, is unchanged. This is consistent with the overall "softness" of the 
approximation being driven mainly by the finite terms in the radiation functions, and not by the choice 
of kinematics map. It is therefore quite possible that one could find another process in which the finite 
terms in the GGG functions would be too large, and the ones in the Ariadne ones just right. 

Returning to fig. 10, as the number of emissions grows, there remains a peak near ME/PS = 1 
(note in particular the logarithmic y axis), but the width of the distribution grows progressively larger, 
indicating that there is a larger number of individual phase-space points in which the pure shower 
is not in agreement with the matrix element. These could be due either to subleading logs, or to 
finite terms in the higher-body matrix elements, not captured by the pure shower. Recall, however, 
that, since we are looking at a flat phase-space scan, we are biased towards the regions of phase 
space where matrix element corrections are important, with the strongly ordered regions occupying 
a progressively smaller volume of the total sampled space. Thus, we interpret this broadening of the 
weight distribution not so much as a sign of any breakdown in the shower approximation itself, but 
rather as illustrating why it is desirable to match to fixed-order matrix elements, a point to which we 
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return below. 



4 Improved Showering 

4.1 Improving the Logarithmic Accuracy: 2 — 3 

A complete set of second-order 2 — 3 (one-loop) and 2 — )• 4 (tree-level) dipole-antenna splitting 
functions is given in ref. [36]. Ultimately, a matching to these functions to all orders in the shower 
would be required to reach formal precision beyond LL, but we note that a simpler, partial matching 
can already be carried out at the 2 — 3 level, to the terms generated purely by the running of the 
coupling. 

These terms can be uniquely identified both in the shower expansion and in the fixed-order antenna 
functions by the fact that they are proportional to the QCD one-loop-running coefficient bo (57), which 
appears in the expansion 




+ 0\ 



a. 



+ 



(87) 



By matching the terms arising from the expansion of as in the shower to the actual feo-dependent 
pieces of the one-loop antenna functions, we obtain a set of universal corrections at the one-loop level 
which stabilize the scale dependence of the resulting calculation to next-to-leading logarithmic (NLL) 
accuracy. 

In particular, we may extract the relevant terms of the one-loop antenna functions by isolating their 
jij-dependent pieces, which are generated purely by quark loops. For gluon emission, the one-loop 
antennae in ref. [36] all contain the following nj-dependent logarithms, 

1 1 ( \ 

rif - (In Vij + In y^k) a^jk " 6 I J ^^'^'^ ' ^^^^ 

where a^^f, denotes the corresponding tree-level antenna function, and p± is defined exactly as in 
Ariadne and Vinci a, i.e., 

pI = . (89) 

SIK 

Because the default renormalization scale used in ref. [36] is 

/^GGG = SIK , (90) 

a redefinition of the renormalization scale from sjk to would absorb the entire term, eq. (88), into 
the definition of the coupling at tree level. Simultaneously, the A'^c'-dependent logarithms generated by 
the same choice can easily be verified to cancel equivalent pieces in the one-loop function, however 
the latter cancellation is not exact due to the presence of additional terms in the one-loop function 
which do not originate from renormalization. (To absorb also these terms would require full one-loop 
matching.) We note that this is nothing but a "renormalization-group-improved" effective redefinition 
of the tree-level coupling which has been known for a long time [53] and is the reason why the default 
renormalization scale both in ViNCiA and in virtually all other Monte Caiios^ is p^. 



''Note, however, that eq. (89) is the only definition of p± that exactly matches the actual 6o-dependent one-loop terms. 
Parton shower models not based on the dipole-antenna picture, which make approximations to this p± definition, will 
therefore necessarily have small bo-dependent remainders left uncanceled. 
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However, rather than just hardcoding one particular choice, we shall here instead interpret it in 
the context of a second-order matching condition, which will allow us the flexibility to estimate the 
remaining uncertainties by varying the scale freely and partly canceling the dependence on it via the 
matching condition. This effectively pushes the effects of the scale variation one order higher in QCD. 

The relevant (partial) matching equation for gluon emission is 

«s(pi)a° fc = a.(/ii,s)P,^^^"aO., , (91) 

which, expanded to first order, easily yields the following form for the scale-stabilizing multiplicative 
factor: 

gluon emission : P^^^^^ = + a^6o In ^ , (92) 

where /ips is the renormalization scale used by the shower evolution and the renormalization scale 
of as in the correction term constitutes an ambiguity of yet higher order. In order to be conservative, 
we wish to make the effects of the scale cancellation produced by this term as small as possible. We 
therefore evaluate the as in eq. (92) at the largest scale in the 2 — > 3 splitting, sjk- Any further 
optimization would amount to a beyond-NLL effect. 

Qualitatively, the scale stabilization works as follows. If fipa, is chosen large, then the correc- 
tion factor, eq. (92), becomes greater than one, hence partly compensating for the lower value. 
Conversely, if a very large as is used at the LL level, the logarithm in the connection term becomes 
negative, and again acts to stabilize the result. We note that, in extreme cases, the correction term 
could in fact become larger than unity. As this would imply a divergent perturbative expansion any- 
way, ViNCiA therefore restricts the range of allowed values to < P^^^m < 2. 

For gluon splitting to quarks, the one-loop antenna functions do not contain universal logarithms 
inp±. Instead, the universal n/-dependent terms are [36] 

n/-ln(y,,-) = n;-ln(-^j , (93) 

where m^g is the invariant mass of the quark-antiquark pair produced in the splitting. In this case, 
we see that the "optimal choice" of renormalization scale is not but m^g. The coiTcsponding 
scale-stabilizing tenii for gluon splitting is therefore 



gluon splitting : Pf^^^' = ( 1 + a.feo In ( ^ ) ) . (94) 



2 



Again, one can easily verify that the iVc -dependent logarithms generated by this choice do partly 
cancel similar pieces in the same one-loop antenna functions, and again, the latter cancellation is not 
exact due to additional pieces unrelated to renormalization. Finally, we note that since 

m\q > pI (95) 

over all of phase space, the effect of this stabilization is to reduce the total number of gluon splittings 
slightly, as compared to what would be obtained without the stabilization terms and using as(pj^) for 
both gluon emissions and gluon splittings, as is traditionally done in shower Monte Caiios. 

These scale stabilizing terms have been implemented in the ViNCiA code for all 2 — )• 3 splittings 
since version 1.020, with an option to switch them on and off to investigate their effects. The default 
in the code is to leave them on. 
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Figure 12: The value of {R4) differentially over 4-parton phase space, with p± ratios characterizing 
the first and second emissions on the x and y axes, respectively. Strong ordering in p± (left) compared 
to no ordering (right). Gluon emission only. Matrix-element weights from MadGraph [51,52], 
leading color (no sum over color permutations). 



Let us emphasize again that this is not a complete one-loop matching. With the scale variation, 
we seek only to evaluate — the scale variation. We do not make any assumption that this variation 
is representative of the entire remaining uncertainty, on which we have several other, independent, 
handles, to which we shall return below. The procedure of employing scale variation alone as a (poor 
man's) estimate of the full uncertainty is obsolete in this framework. 

4.2 Improving the Logarithmic Accuracy: 2 — 4 

While parton emission using trial branchings can easily be made to cover the full phase space for a 
single emission, the same is not true for multiple emissions. Due to the requirement of strong ordering, 
some regions of phase space may be inaccessible, leading to so-called dead zones. This also happens 
in strongly ordered dipole-antenna showers, for example in regions where several emissions happen 
at closely similar emission scales, as shown in ref. [34, 54]. Since gluon emission and gluon splitting 
processes have different singularity structures and are treated slightly differently, we first consider 
the dominant case, that of gluon emission. We then give a few brief remaifc about gluon splitting, 
although we defer most of the details of that discussion to another publication [55]. 

4.2.1 Gluon Emission 

A plot from ref. [34], showing the dead zone for Z — > qggq in a pj^-ordered dipole-antenna shower, 
is reproduced in the left-hand pane of Fig. 12. Each bin of this 2D histogram shows the value of i?4, 
eq. (83), averaged over all 4-parton phase space points that populate that bin. The black zone above 
the strong-ordering line corresponds exactly to the spike on the left-hand edge of the plots in fig. 10 
(the underflow bin). 

If one simply removes the strong-ordering condition, equivalent to ordering the emissions only by 
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the nesting of the factorized phase spaces, the dead zone obviously disappears. However, this comes 
at the price of introducing a substantial double counting, which also extends deep into the ordered 
region. We illustrate this on the plot in the right-hand pane of fig. 12, in which {R4) S> 1 both in the 
unordered region as well as in parts of the ordered region where the agreement was previously good. 
Clearly, therefore, the ordered description, in the left-hand pane, is a better overall approximation to 
QCD, even if it does include a dead zone. 

In order to go further, we must understand the physics behind the ordered and unordered approx- 
imations to the matrix elements. Why is ordering so important? The first exploration of this goes 
back to the early eighties. Then, it was realized that parton showers ordered only in parton virtuality 
(which is equivalent to a pure phase-space ordering in that language) represent an essentially incoher- 
ent addition of independently radiating monopoles. In phase-space regions where the contributions 
from each such monopole term are comparable, interference effects can become large. Without them, 
the pure phase-space ordered shower gives a substantial overcounting of those regions, as compared 
to matrix elements [28,29]. As Marchesini and Webber showed [40], this double counting can be 
approximately identified with temis conesponding to non-angular-ordered emissions, and hence the 
procedure to impose coherence on traditional parton showers has since been to impose such an order- 
ing, either implicitly by the choice of evolution variable, as in Herwig, or explicitly as a veto on the 
generated trial emissions, as in Pythia. 

In dipole-based shower models, soft coherence inside each dipole is guaranteed, regardless of 
the ordering variable, by using dipole-based radiation functions instead of the DGLAP ones, but the 
problem still exists; it has just been pushed one order higher in the number of interfering partons (see, 
e.g., refs. [34,49]). With pure phase-space ordering, dipole-antenna showers essentially represent an 
incoherent addition of independently radiating dipoles. An independent addition of two such dipoles 
would result in a substantial overcounting in all regions where several such dipole terms contribute 
simultaneously at similar levels, i.e., in regions where dipole-dipole interference effects (or, equiva- 
lently, multipole effects) would be important. Again, it would be interesting to investigate whether 
some variant of angular ordering could be used to restore a more coherent behavior, but in the context 
of the dipole-antenna formalism we develop here, we have not been able to find such a solution. In 
part, this owes to a strict ordering in angle having some disadvantages in our language; being frame- 
dependent, it would not respect the Lorentz-invariant dipole phase-space factorization we rely on, and 
it also classifies a subset of infinitely soft and/or collinear emissions as happening at finite values of 
the evolution variable, which would lead to ill-defined evolution integrals, see, e.g., ref. [34]. 

Instead, let us recall the basic motivation for angular ordering: to approximately remove the dou- 
ble counting caused by incoherent addition of interfering diagrams of similar magnitudes. In the 
parton shower language, each such term is associated with a divergence in the energy times angle 
of the emission. In the region where several terms would nominally be large, the angular ordering 
requirement forces at most one of them to contribute — approximating the destructive-interference 
effects by killing the non-ordered contributions. In dipole-based approaches, however, the leading 
divergence of the gluon radiation functions occurs unambiguously in the of the emitted gluon. It 
therefore seems sensible to use p_\_ as the measure for the ordering, and thereby implicitly for the size 
of each of the contributing terms. As can be seen from the left-hand pane of fig. 12, an ordering in 
p_L, yields a quite good average approximation to the full 2 — 4 matrix elements over most of phase 
space. 

We included the above discussion to motivate that, while there is an important physics issue behind 
strong ordering and also behind the choice of functional form of the ordering variable, there is nothing 
particularly important about imposing it as a step function in that variable. On the contrary, the actual 
destructive-interference terms in the matrix elements exhibit a smooth behavior. To remove dead 
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Figure 13: The value of {R4) differentially over 4-parton phase space, with pj_ ratios characteriz- 
ing the first and second emissions on the x and y axes, respectively. Smooth ordering in p± (left) 
compared to smooth ordering in m£) (right). Gluon emission only. Matrix-element weights from 
MadGraph [51,52], leading color (no sum over color permutations). 



zones to all orders while simultaneously improving the shower approximation also in the ordered 
region, we therefore propose to change the condition of strong ordering to a smoother condition with 
the same limiting behaviors. 

Specifically, while we retain strong ordering as an option in ViNCiA, by default we replace the 
strong-ordering condition of conventional parton showers in gluon emission by the smooth suppres- 
sion factor ^ 

Gluon Emission : Bord ^ PimnP^^ = 9 P^^ > (96) 

Pi +pi 

where p± is the smallest p± scale among all the color-connected parton triplets in the parent config- 
uration (i.e., a global measure of the "current" p± scale of that topology); p'i is the scale of the trial 
2 — )• 3 emission under consideration; and P^^ is defined in eq. (82). 

Since the antenna function for the previous branching is proportional to 1 /p^ , the net effect of 
this term, in the unordered region, is to replace that divergence by a damped factor, l/{p]_ + p^)- 
The correction is thus constructed such that P^^ remains unmodified in the strongly ordered limit 
p± p±, and therefore will not affect the leading-logaiithmic behavior of the parton shower. It then 
drops off to ^P^^ for p± = p±, and finally tends smoothly to zero in the limit of extreme unordering, 
p± > p±. 

The ratio of the resulting shower to matrix elements is shown in the left-hand pane of fig. 13. 
Comparing this distribution with those in fig. 12, we indeed see that not only has the dead zone been 
removed, without introducing any serious overcounting of it, but the quality of the approximation has 
also been improved inside the ordered region. 

For completeness, in the right-hand pane of fig. (13), we also show how the approximation would 
have looked if the alternative measure m|, = 2 min (m?j ,m|^) had been used instead of p± in the 
suppression factor eq. (96). Although there is still clearly an improvement over the pure phase-space- 
ordered case — the dead zone has been eliminated — it is much less convincing than for p±, as the 
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Figure 14: Smoothly ordered parton showers compared to matrix elements. Distribution of 
log]^Q(PS/ME) in a flat phase-space scan. Contents normalized by the number of generated points. 
Gluon emission only. Matrix-element weights from MadGraph [51,52], leading color (no sum over 
color permutations). Compare to fig. 10 for strong ordering. 



weights are larger in the region above the thin horizontal red line. 

To illustrate how this approximation evolves with parton multiplicity, we show the distribution of 
the log of the PS/ME ratio with this modification, in fig. 14, for Z — )• 4, 5, and 6 partons, including 
only leading-color gluon emission. One observes a marked improvement with respect to the strongly 
ordered approximations, fig. 10, for all multiplicities. In particular, not only the dead zones but 
also the large tails towards low PS/ME ratios visible in the higher-multiplicity plots in fig. 10 have 
disappeared, which we interpret as a confimiation that the logarithmic accuracy of the approximation 
has indeed been improved. Notice, however, that the Ariadne functions (where we have here used 
the TpAR kinematics map for both qg and gg antennae, hence the explicit label on the plot) still tend to 
shift the shower approximations systematically towards softer values, whereas the GGG ones remain 
closer to the matrix elements. 



4.2.2 Gluon Splitting 

For gluon splitting, there is no soft singularity, only a collinear one. This means there is now only a 
single log-enhancement (instead of a double log) driving the approximation and competing with the 
(uncontrolled) finite terms. It is therefore to be expected that the LL approximation to gluon splitting 
is significantly worse, over more of phase space, than is the case for gluon emission. 

Furthermore, if the two neighboring dipole-antennse that share the splitting gluon are very un- 
equal in size, e.g., as a result of a preceding close-to-coUinear branching, then higher-order matrix 
elements and splitting functions unambiguously indicate that the total gluon splitting probability is 
significantly suppressed. This is not taken into account when treating the two antennae as independent 
radiators. This effect was already noted by the authors of Ariadne, and a first attempt at including it 
approximately was made by applying the following additional factor to gluon splittings in Ariadne, 
in addition to the strong-ordering condition, 

Gluon Splitting (Ariadne) : Oord ^ QordPAviP^^ = ©ord — — — P^^ , (97) 

SIK + SN 

where is the invariant mass squai^ed of the neighboring dipole-antenna that shares the gluon 
splitting, and sjk is the invariant mass squared of the dipole-antenna in which the splitting occurs. 
The additional factor reduces to unity when the two neighboring invariants are similar; it suppresses 
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splittings in an antenna whose neighbor has a very small invariant mass; and it slightly enhances 
splittings in dipole-antennae whose neighbors have very large invariant masses. This modification is 
ad hoc, but formally is beyond LL, and thus does not spoil the shower's properties at this order. In 
practice, it greatly improves the shower's approximation to matrix elements. 

For our purposes here, we first replace the O function in Ariadne by the same smooth unordering 
suppression factor we used above, 

Gluon Splitting : OordPlL ^ PimpPAnPLL , (98) 

which gives an overall approximation at least as good as that of Ariadne, without any dead zones. 
The overall agreement between even this improved gluon splitting and explicit matrix elements is 
still far from perfect, however, due to the intrinsically smaller relative size of the logs driving the 
approximation. We are in the process of preparing a more dedicated study of this issue, including 
quark mass effects [55], and hence defer further discussion of this topic for the time being. For later 
convenience, we define this adjustment factor to be unity for gluon emission, 

Gluon Emission : Pau = 1 , (99) 
Gluon Splitting : ^Ari = — — — • (100) 

SIK + SN 



4.3 Subleading Color 

In the dipole-antenna formalism, a general result [56] is that the subleading-colour effects in a single 
qgg---gq chain can be taken into account by including a subleading antenna spanned directly between 
the q and q associated with a color factor —l/N^ relative to the leading-color antennae (which are 
proportional to Ca)- 

However, in the context of a probabilistic framework, such as shower Monte Carlo algorithms, the 
negative sign of this antenna means it cannot be treated on the same footing as the (positive-definite) 
LC ones^. Moreover, it is not possible to define a unique LC color assignment to the emissions 
generated by it, and hence the subsequent shower evolution (and the infinite-order approximations 
generated by it) would be ill-defined. 

Instead, in traditional parton-shower applications, this correction is partly treated by associating 
quark emission terms with Cp instead of Ca, thereby correctly absorbing the coUinear singularities 
of the correction term into the LC ones, at the price of introducing a subleading-color ambiguity in 
the soft singularity structure. To improve on this, one could imagine, e.g., trying to be more clever 
about in which phase-space regions to use Cp and in which Ca even for emissions off gluons [57]. 
But in both cases the simplicity of the correction term would then be less explicit. 

Instead, we here attempt to reabsorb the subleading-color correction systematically into the leading- 
color antennse, using a smooth partitioning, which integrates to reproduce the double poles of the 
corresponding subleading-color one-loop antenna functions. 

Consider the evolution integral off an n-paiton configuration in the massless approximation: 

/ dsij dsjk 2 ' (101) 

IK&LC-' ^IK 

'One could in principle imagine flipping the sign of the event weight when generating emissions with it, but such a 
procedure would drive the convergence rate of the resulting algorithm to become infinitely slow at asymptotically large 
energies. 
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where a is a dimensionless antenna function, and the scaled invariants are defined by yij = Sij/sijk = 
Sij/sjK, and where the sum is over all color-connected pairs, that is all LC antennse. Changing 
integration variables to dimensionless quantities, the integration measure becomes independent of the 
size of the antenna phase spaces and we may replace the sum of integrals by the integral of a sum, to 
which we can add a subleading-color piece, 



Jdyijdyjk ^iK{yij,yjk) - ^aNLc{yij,yjk) 



.IKeLC c 



(102) 



Finally, we may partition the coiTcction term among the LC pieces by introducing a partitioning 
function, /, 



E~ , X flK aNLc{yij,yjk) 
aiK{yij,yjk) 1 - ^ '—^ 



IKeLC 



Nc aiK{yij,yjk) 



Pnlc 



(103) 



where the underbraced tenii can now be implemented straightforwardly as a veto probability, Pnlc^ 
in the shower evolution. 

For / to be a consistent partitioning, it must give unity when summed over all the LC antennae. 
The prescription we use is to absorb corrections into each term in proportion to the relative size of 
that terni, i.e., 

f _ aiKiyij,yjk) 

JIK - ^ : 7 r ■ (104) 

l^ABeLC aAB{yar,yrb) 

For the functional form of a, we give two options, 

1. Only the Eikonal part of the qq antenna function is included in a j\fLc^ corresponding to 

/ N '^{'^ - Vij - Vjk) 

aNLc[yij,yjk) = — (105) 

Vijyjk 

This integrates to give the correct poles of the l/N^, piece of the one-loop antenna func- 
tion (called Al in ref. [36]). 

2. The full qq antenna function (called in [36]) is included in clnlc^ corresponding to 

( \ '^{'^ - yij - yjk) , yij , Vjk 
aNLc{yij,yjk) = — + — + — (106) 

Vijyjk Vjk yij 

This integrates to give the correct 1/e^ and l/e poles of the l/N^, piece of the one-loop antenna 
function (called ^3 in ref. [36]). 

For the corrected emission probability to be positive definite, the condition 

aNLC < -j^aiK (107) 

JIK 

must be fulfilled. The Eikonal terms are guaranteed to respect this by a wide margin, but there 
can in principle be subleading differences between a//< and a^Lc at the level of single poles and/or 
finite temis, both of which can be influenced to some extent by user-controlled settings in ViNClA. 
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Figure 15: LC and NLC parton showers compared to matrix elements. Distribution of log^o (PS/ME) 
in a flat phase-space scan. Contents normalized by the number of generated points. Gluon emission 
only. Matrix-element weights from MadGraph [51, 52], full color (summed over color permuta- 
tions). Compare to figs. 10 and 14. Note: the black histogram here does not look identical to the one 
in fig. 10; this is because the distribution shown here is after summation over all color-permutations. 
Also, the shower expansion here uses Ca for both qg and gg emission antennse. The same is true for 
the thin solid line, as compared to the solid histogram in fig. 14. 



Although we have not encountered any problems with such corrections becoming large in our own 
studies so far, we have inserted a numerical safeguai^d in the code, limiting the size of each NLC 
correction to be at most half of the corresponding LC term. 

We compare the shower expansion corrected as in eq. (103), using the definition of flTVLC given by 
eq. (106), to color-summed matrix elements in Fig. 15, for Z — )• 5 and Z — 6 partons (gluon emission 
only), for both the strong- and smooth-ordering options. For reference, we show the pure LC shower 
expansions as well, as in figs. 10 and 14, although we here set both the qg and gg color factors equal 
to Ca, to better illustrate what happens when only the NLC correction is switched on and off. Note 
also that the previous comparisons in this paper were made to leading-color matrix elements for each 
individual color structure separately. Such a comparison is not meaningful once color interference 
effects ai^e taken into account, and hence we are here using the full color-summed matrix elements, 
and are likewise summing over permutations of the gluon momenta in the shower expansion. 

Firstly, we notice that the LC shower expansions do indeed overcount the matrix elements if we 
just use Ca for both qg and gg antennse, as expected (the solid filled histogram and thin solid Unes 
are shifted to the right of the vertical dot-dashed line that represents perfect agreement). When we 
switch on the NLC coiTcction in the manner described above, both the strongly and smoothly ordered 
approximations aie noticeably improved, cf. the blue dashed and thick solid cui"ve, respectively. The 
NLC correction is therefore switched on by default in the ViNCiA code. 

The procedure described above is correct for the case of a single qgg---gq chain and then correctly 
takes into account the infrared singularities arising from the first subleading-color term, which is 
proportional to l/N^;. We refer to this modification as "next-to-leading color" or NLC. (Strictly 
speaking, there are also corrections proportional to l/Nc, but these are akeady taken into account 
in leading-color showers by including g ^ qq splittings.) However, due to the inherent ambiguity 
in assigning a color flow to these corrections even in the soft limit, we cannot be certain that further 



38 



subleading terms are also correctly described. I.e., we do not expect to reproduce the coiTcct 
terms in all singular limits. 

When there are several colour chains, such as after a g ^ qq splitting (and hence already sup- 
pressed by at least 1/Nc), we generalize the treatment to each chain separately. That is, we do not at 
this point attempt to include interference effects between different color-singlet systems. 



5 Matching 

In this section, we describe our strategy for incorporating a detailed matching to tree-level matrix 
elements. The philosophy is similar to that pioneered by Sjostrand in refs. [28,29], and hence also to 
the POWHEG formalism, but we here generalize the method to include arbitrary-multiplicity tree-level 
matrix elements. The inclusion of the NLO virtual corrections to the lowest multiplicity was treated 
in ref. [33] for an arbitrary tree-level matching strategy. 



5.1 Matching Strategies 

Given a parton shower and a matrix-element generator, there are fundamentally three different ways 
in which we can consider matching the two: 

1. Unitarity: The oldest approach [28,29] consists of working out the shower approximation to a 
given fixed order, and correcting the shower splitting functions at that order by a multiplicative 
factor given by the ratio of the matrix element to the shower approximation, phase-space point 
by phase-space point. We may sketch this as 

Matched = Approximate ^^^^^ . (108) 

Approximate 

That is, the shower approximation is essentially used as a pre-weighted (stratified) all-orders 
phase-space generator, on which a more exact answer can subsequently be imprinted order 
by order in perturbation theory. In our notation [33], this translates to applying the following 
coiTcction factor to each antenna function Cj (or any other kind of shower splitting kernel) 

where the sum over j runs over all possible ways the shower could have generated the n- 
parton state from n — I partons^. So long as the adjustment factors are less than or 

equal to one, they can be interpreted as probabilities, and the adjustment can be accomplished 
by means of the veto algorithm Monte-Carlo technique. This constraint can essentially al- 
ways be satisfied through appropriate choice of the finite terms in the antenna functions Oj. 
When these coiTcction factors are inserted back into the shower evolution, they guarantee that 
the shower evolution off n — 1 partons con^ectly reproduces the n-parton matrix elements, with- 
out the need to generate any separate n-parton samples. Moreover, since the corrections modify 
the actual shower evolution kernels, the corrections are resummed in the Sudakov exponen- 
tial, and finally, since the shower is unitary, an initially unweighted sample of (n — 1) -parton 



Note, however, that this gets substantially more complicated if the shower process is not completely Markovian, a point 
we shall return to. 



39 



configurations remains unweigiited, with no need for a sepai^ate event-unweigliting or event- 
rejection step. (Technically, the exponentiation allows beyond-LL corrections to be resummed, 
thus improving the logarithmic accuracy of the result, while the explicit constraint of unitarity 
ensures that the additional non-logarithmic terms that are also exponentiated by this procedure 
do not lead to disasters.) There are thus several quite desirable features to this kind of match- 
ing strategy, which is cuiTcntly employed by Pythia, Powheg, and ViNClA. However, since 
traditional shower expansions quickly get more complicated as a function of the number of 
emissions, this strategy had only been worked out for a single additional emission prior to this 
paper (although the Menlops strategy [26] does allow to combine a unitary matching of the 
first emission with traditional non-unitary methods for multi-jet matching). Below, we shall 
generalize the unitarity method to arbitrary multiplicities and, as a proof of concept, present 
a concrete implementation spanning four successive emissions, including all subleading color 
terms. 

2. Subtraction: Another way of matching two calculations is by subtracting one from the other 
and correcting by the difference, schematically 

Matched = Approximate + (Exact — Approximate) . (HO) 

This looks very much like the structure of an NLO fixed-order calculation, in which the shower 
approximation plays the role of subtraction terms, and indeed this is what is used in strategies 
like Mc@NLO [22, 58, 59]. In particular since eq. (110) appears much simpler to the fixed- 
order community than eq. (108), this type of approach has received much more attention than 
the unitarity-based one above (though, to be fair, the PowHEG [30] approach represents a kind 
of hybrid between the two). In this approach, the corrections are not resummed; the events are 
not unweighted — we can even have negative weights, at phase-space points where the approx- 
imation is lai^ger than the exact answer; and we need a separate phase-space generator for the 
n-parton coixection events. And finally, like for the unitarity-based case above, a systematic 
way of extending this strategy beyond the first additional emission was not previously avail- 
able. All these issues are, however, less severe than in ordinary NLO approaches, and hence 
they are not viewed as disadvantages if the point of reference is an NLO computation. Since 
the correction terms are applied by adding (or subtracting, depending on the sign of the weight) 
events, we refer to this type of matching strategy as subtraction. 

3. Slicing: The last matching type is based on separating phase space into two regions, one of 
which is supposed to be mainly described by hard matrix elements and the other of which is 
supposed to be described by the shower. Basically, this amounts to a subtractive approach 
in which the shower approximation is set to zero above some scale (effectively a dead zone 
is forced on the shower by vetoing any emissions above a certain matching scale), causing 
the matched result to be the unsubtracted matrix element in that region, modulo higher-order 
collections. 

Matched (above matching scale) ~ Exact (1 + 0{as)) , (HI) 

and since the leading behavior of the matrix elements and the shower approximation are as- 
sumed to be the same below the matching scale, the small difference between them can be 
dropped, yielding the pure shower answer in that region. 

Matched (below matching scale) = Approximate -|- (Exact — Approximate) ~ Approximate . 

(112) 
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Since this strategy is discontinuous across phase space, a main point here is to ensure that the 
behavior across the matching scale be as smooth as possible. CKKW showed [6] that it is possi- 
ble to remove any dependence on the matching scale through NLL precision by careful choices 
of all ingredients in the matching; technical details of the implementation (affecting the 0{as) 
terms in eq. (Ill)) are important [13]. The MLM [4, 18] approach is also an example of this 
type of matching. However, we note that the slicing strategy inherits almost all of the problems 
that pure subtraction has; the corrections are not resummed, the events are not unweighted (but 
at least we avoid the negative-weight issue), and we need a separate phase-space generator 
for the n-parton correction events. In addition, the dependence on the unphysical matching 
scale has appeared, which in general is non-vanishing and may be larger than NLL unless the 
implementation matches the theoretical algorithm precisely [13]. However, due to the work of 
CKKW and others [6, 13, 14, 18], a systematic way of extending this strategy beyond the first 
additional emission is available, and a strategy has even been proposed whereby this matching 
type could be extended to NLO precision [60]. The Menlops approach [26] is also available 
to combine it with POWHEG. Since this type is already well developed, therefore, we shall not 
consider it further in this paper, but note that one could still obtain it from our formulae for 
subtractive matching (eq. (22) of ref. [33]), by inserting the appropriate phase-space cutoffs at 
the matching scale. For this reason, we refer to this strategy as scale-based or slicing. 

To summarize, in this paper, we focus on the extension of the unitary matching strategy to arbitrary 
numbers of emissions at tree level. We shall also include an NLO matching to the Bom multiplicity 
using the prescription from ref. [33]. 

5.2 Matching to Tree-Level Matrix Elements: Leading Color 

The formalism we shall describe here represents a generalization of the one presented in ref. [33]. It is 
based on using the trial generator described in the previous section as a phase-space generator, whose 
phase-space weights can be expanded to the required order and compared to the matrix-element an- 
swer at the same order. A correction can then be imposed before generating the next trial emission. 
For comparison, most other approaches are cunently based on generating separate event samples for 
different jet multiplicities and then post facto attempting to remove the overlaps ("double counting") 
between them. Here, we generate one sample ab initio, where every event starts at the lowest multi- 
plicity and is then successively matched up to the desired orders. 

Compared to shower evolution, matrix-element (ME) evaluations are computationally intensive. 
It is therefore desirable that the matching algorithm involve as few ME evaluations as possible. In an 
ordinary shower approach, the effective weight of each phase-space point depends on all the possible 
shower histories that could contribute to it, resulting in a factorially increasing dependence on the 
multiplicity (for each color configuration). In CKKW-type approaches, this is partly circumvented by 
always selecting only one history, the "most singular one" according to the kt algorithm. Since this 
does not exactly correspond to an inversion of the parton shower, the matching only really addresses 
the LL overlaps, and the higher-order discrepancies ai^e "removed" by introducing an explicit cut on 
the phase-space region in which matching is applied, by the so-called "matching scale", which limits 
the numerical size that any subleading divergence could attain. 

In contrast, since we match at each successive order, each emission (up to the matched orders) will 
necessitate at least one matrix-element evaluation as well as the corresponding shower weight, which 
in turn will involve matrix-element evaluations at the preceding orders. At first sight, this may sound 
extremely expensive. Note, however, that we are doing the matching of "all the samples" once and for 
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all, so that only one "run" will be necessary, rather than a separate one for each multiplicity. (Separate 
"runs" for different multiplicities would then spread a comparable number of matrix-element evalu- 
ations across the different "runs".) In that context, the scheme should not be more expensive than 
current ones, provided that a formalism can be found that minimizes the number of matrix-element 
evaluations that are still necessary to determine the trial weight. 

One sufficient condition for minimizing the number of required matrix-element evaluations is the 
Markov condition: each shower step should depend only on the current n-parton configuration and 
not on its previous history. This in turn implies that, in order to compute the trial weight, only the 
histories one step back have to be considered, rather than all possible clusterings all the way to the 
Bom. As mentioned above, this would not be true of ordinary strongly-ordered parton showers, where 
the restart scale for each configuration would depend on which parton was the last emitted one. 

A simple prescription that does obey the Markov condition is to generate trial emissions for every 
antenna in the n-parton configuration over their full phase space, irrespective of the current ordering 
scale. Without matching, this would lead to a large overcounting, as was illustrated in fig. 12, but with 
matching, the total shower weight can be calculated and the coiTcsponding matrix-element coiTcction 
made, with two added benefits: 1) the removal of the strong-ordering condition explicitly prevents 
any dead zones from appearing in the trial space, and 2) since the trial weight generated this way 
will represent an overestimate, it will be possible to impose the matching by multiplying by a factor 
smaller than unity, which can be translated into a probabilistic veto of the trial branching. 

Expanded to tree level (all Sudakovs set equal to unity, fixed a^), the trial generator will produce 
the following total weight for a specific color-ordered point in (n + 1) -parton phase space when 
summed over possible contributing n-parton ones, 

«;J;ti(M.+l) = $:«trial-,({Pn}[^] ^ {rfn+l) | M^J^ ({p} ^1 ) | ^ , (113) 
j 

where {p}n+i represents the color-ordered momenta of the (n + l)-parton state, j mns over the 
possible n + 1 — )• n clusterings, and {p}^n represents the color-ordered set of n momenta obtained 
by the j'th 3 — 2 clustering of the (n + l)-parton state according to the selected kinematics map. It 
is important to note that | M^^^^ \ ^ here represents the tree-level squared amplitude for the particular 
color configuration under consideration, i.e., without any color averaging performed. 

The improved matching to smoothly ordered LL antenna functions described in sec. 4.2.2 merely 
consisted of multiplying otriai-j in eq. (113) by the LL smooth-ordering accept probability, thus 
replacing the trial factors by their LL counterparts, 

LL Matching : atriai-j ^ Oj^ = P^^ept-j «triai-i , (114) 

where P^^^pt-j expresses the unmatched shower accept probability, including the LL acceptance 
probability, P^^, of eq. (82) and the improvement factors Pimp of eq. (96) and Pau of eq. (97), 

-^accept-j ~ PimpPAriP • (115) 



We shall now apply a final multiplicative accept probability, P;^^, defined such that it takes us 
from the approximation that would have been generated by P^^ept-j ^lone to the full matrix element. 
It has the simple definition 

pME ^ |M.+i({pWi)P 

E.aP(Ml'V{pWi)|M„({p}W)P 
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where a^^ was defined in eq. (114). By summing over all shower histories, this can easily be verified 
to generate the correct total weight. Note also that, the ME accept probability does not have any 
dependence on j and is thus the same for all contributing n — )■ (n + 1) branchings. 

Note that i^^^ involves one evaluation of |M„_|_i({p}„_|_i)p and one evaluation of each of the 

reclustered configurations, |M„({p}n')p (one for each possibility for [k], of which one has already 
been evaluated as part of the current matched shower history), and hence the total number of matrix- 
element evaluations required at each order grows linearly with the multiplicity, rather than factorially 
as would have been the case for a non-Markov evolution. 

Note also that, since we make no distinction between "shower events" and "correction events", we 
may use the shower Os as a common prefactor on all the accept probabilities. There is thus only one 
as for each history. Due to the scale-canceling partial one-loop matching discussed in sec. 4.2.2, the 
default, to one-loop accuracy, is thus to use as{p_L) for all the terms, regardless of other scale choices 
made in the generator. 

Another convenient way of writing P^^i is the following. 



pME -j^ _j_ I — /i.-r-L \ LJ- J it-r-L / I ~K VLJ-j'i- Li^ji'-r±/\ — ii-vLJ^j" /I (117) 

where the connection to the subtraction approach described in ref. [33] (as well as to other subtrac- 
tion schemes) becomes readily apparent, since the numerator in eq. (117) is nothing but a shower- 
subtracted matrix element. 

In order to transform the unitary strategy described here to a non-unitary subtractive one, it would 
therefore suffice to apply the factor, eq. (117), as an event reweighting, rather than as a branch-accept 
probability. The events then do not have unit weights any more, and a subsequent unweighting step 
would be necessary, as in other subtractive approaches. The fact that we are here doing the matching 
phase-space point by phase-space point, however, means that we here have the in front, which 
should mean that even the subtraction-based version could never generate negative weights. In current 
subtraction approaches such as Mc@NLO, the is generated as a separate sample, and in that case, 
the collection term by itself can of course yield negative collection events. The cancellation here is 
more elegant and not only yields positive- weight events but is also better protected from fluctuations; 
in Mc@NLO and POWHEG, each event sample is uncorrelated and therefore the phase-space point of 
an event will never be hit exactly by the counter-events. In limited-statistics samples, it is therefore 
events with slightly different momenta that have to compensate each other, whereas the proposal here 
achieves an exact cancellation in one and the same phase-space point, event by event. Thus, instead 
of having one event with weight -i-3 and one with weight - 1 , we would here simply get one event with 
a total weight of +2. 

A final technical note is that the MadGraph matrix elements must be evaluated on-shell, and 
hence one must first set the value of the MadGraph Z mass equal to ^/s, even if this is not equal to 
the physical Z mass; the important thing is that the incoming momentum be on-shell. We also set the 
value of the strong coupling gs = lin MadGraph, equivalent to factoring it out of the problem. 



M„+i({pWi)P -Efc«r'({p}l" ^ {pWi)|M„({p}W)|2 



5.3 Matching to Tree-Level Matrix Elements: Subleading Color 

When summing over all events, the full answer contains an averaging over all permutations of color 
orderings in every phase-space point. In event generators, two different color structures in one and 
the same phase-space point are viewed as two different events with the same color structure, but in 
two different phase-space points. Via the matching above, this sum now reproduces the leading-color 
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tree-level color-averaged matrix element squai^ed. We wish to extend this to include the subleading- 
color contributions as well. Obviously, these cannot be associated with any particular color structure, 
and we must therefore here match across events in different phase-space points (or, equivalently, the 
same momentum configuration, but different color orderings). 

When matching to subleading color, we shall use the specific structure of matrix elements gen- 
erated with MadGraph and imported into the ViNCiA code. The color structure of these matrix 
elements is cast as a color matrix whose diagonal entries form the leading-color contributions. Each 
diagonal entry corresponds to one particular color ordering, hence summing over all the diagonal 
terms and dividing by the number of rows (= number of permutations) is equivalent to averaging over 
colors, in the leading-color limit. The leading-color matching described above thus corresponds to 
matching to a MadGraph matrix element with only one diagonal term being non-zero in the entire 
color matrix, representing the particular color stnicture used in the matching. 

In order to include subleading color, a simple and sufficient prescription is to compute also the 
full color-summed matrix element and include a fraction of it in the matching to each color structure, 
by modifying the LC matrix elements in eqs. (116) and (117) as follows, 

|/Vf.|2 

|M,|2 ^ m' + ^ n\r\2 T.M,M; (118) 

= i^TPitrj • 

where Mi is the amplitude for one specific color ordering. The numerator in the last parenthesis is 
just the full color-summed matrix element squared, and the denominator is the corresponding leading- 
color one. Since both of these are well-defined leading-order physical matrix elements, the term in 
parenthesis is positive definite and hence cannot generate negative weights. 

Note that, in ViNClA's interface to MadGraph, we have so far been using the form in eq. (118), 
since this is the fastest in the context of that particular implementation. In order to compare between 
the leading- and full-color cases, we have implemented an option to switch the subleading corrections 
off, although they are on by default in the program. 

In fig. 16, we show the weight ratios discussed earlier (which are essentially just the inverses of 
for Z ^ 5 and Z ^ 6 partons, now including matching at each preceding order. For the 
shower approximations, we use the default smoothly ordered NLC -improved GGG antennse, with 
three different kinematics maps (solid histogram, thin solid line, and dashed lines, respectively). We 
also compare to the same settings as the solid histogram but using the Ariadne radiation functions 
instead of the GGG ones (thick solid lines). Comparing these distributions to those in fig. 14, we 
see that all the shower models reproduce the matrix elements very well, and hence the differences 
between the shower models are largely canceled by the matching to the preceding orders, as expected. 
At each order, now only a relatively well-controlled and stable matching correction remains, which 
does not appear to exhibit any significant deterioration order by order. Note that we have not applied 
any phase-space cuts here, and hence we find no evidence for any remaining subleading divergences 
in the matrix elements. This is in sharp contrast to slicing- or subtraction-based approaches, where a 
non-zero matching scale is obligatory beyond the first matched order. 
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Figure 16: Smoothly ordered matched parton showers (unordered showers using the improvement 
factors Pimp and PAri) compared to matrix elements. Distribution of log^g (PS/ME) in a flat phase- 
space scan. Contents normalized by the number of generated points. Gluon emission only. Matrix- 
element weights from MadGraph [51,52], full color (summed over color permutations). Compare 
to the unmatched shower distributions in figs. 10, 14, and 15. 



5.4 A note on MadGraph and GGG color factor normalizations 

As a final remark, we note that the subleading-color terms are not uniquely defined. Obviously, if 
the leading-color pieces are not normalized the same way in two different approaches, the sublead- 
ing terms must likewise appear different. This, e.g., leads to some apparent differences between 
MadGraph and the GGG antennae. With color and coupling factors, the MadGraph-GGG corre- 
spondence for the Z — qggq antenna is: 

gtA^^^ {0,1,2,3) = 2|^4Lc(0,l,2,3)p 

where the factor 2 on the MadGraph matrix element cancels the color averaging factor which is 
already present in |M4lcP. which represents a MadGraph matrix element with only one element 
non-zero in the color matrix, the one corresponding to the (0, 1,2,3) color flow squared. In paiticulai\ 
note that the LC coefficient in MadGraph comes with Cp, whereas, in order to construct the fuU 
answer, i.e., including subleading color, the GGG antennae should be combined in the following way. 



- gtCpCA i-a^i0,l,2,3) + -aii0,2,l,3) - j^A^{0,l, 2,3) ] . (122) 



A direct comparison between what would be called subleading color by GGG and by MadGraph, 
respectively, would thus yield a different answer, simply because the piece called LC is normalized to 
CaCf in the former and to Cp in the latter. To verify these normalizations, the validity of eqs. (121) 
and (122) was tested numerically on a large number of phase-space points. 
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6 Uncertainty Bands 



A calculation is only as good as the trustworthiness of its uncertainty bands. Traditional methods 
for evaluating shower uncertainties range from simple comparisons between different models to more 
elaborate variations of salient model pai^ameters within some theoretically or phenomenologically 
justified ranges. 

The former kind is, at best, indicative, but can also be grossly misleading. As a classic example, 
consider two different parton showers with a cutoff at some factorization scale. They would both agree 
there are no jets above that scale, even though a matrix-element-based calculation would certainly 
produce jets in that phase-space region. Comparisons of the Herwig — Pythia kind are therefore 
of little value when pursuing rigorous uncertainty estimates. 

Systematic variation of salient model parameters obviously gives a more trustworthy idea of the 
overall uncertainty, and can also give infomiation about which particular sources dominate. However, 
it requires more careful preparation and more expert input to set up: which parameters to vary, within 
what ranges, and how to make sure the variations are done consistently when combining many tools 
in a long chain of event generation. It also requires substantially more time and resources: for each 
variation, a new set of events must be generated, matched, unweighted, and possibly passed through 
detector simulation. Finally, the ability of a single model to span all possible variations is often limited 
— similarly to above, you still cannot use a strongly ordered shower to estimate what the uncertainty 
associated with the strong-ordering condition itself might be. There is also no way that, for instance, 
Pythia's shower model could be varied to obtain an estimate of what an angulai^-ordered shower 
would give. 

Here, we propose to combine the flexibihty of the ViNCiA formalism to take into account different 
ordering variables, radiation functions, etc., with a treatment of uncertainties that only involves the 
generation of a single event sample, with a time requirement that is not greatly increased compared 
to the case without uncertainty variations. We shall also automate the expert input to some extent, 
reducing the number of choices the user must make. 

The key question to ask is: if we use (matched) parton shower model A to generate a set of 
unweighted events, what would the weight of each of those events have been if we had instead used 
parton shower model B to generate them? By answering this question, we can essentially use any 
parton shower model as a "phase space primer", provided it is still reasonably physical and that it 
does not have any dead zones, and then compute alternative weights /or the same events for any other 
set of assumptions. 

The most trivial part is to note that, if a particular shower model uses Osiai as its radiation function 
for a particular branching, the same branching would have happened with the relative probability 

P2 = Pi , (123) 

in a different model that uses as2 as its coupling (e.g., with a different renormalization scale or 
scheme) and 02 as its radiation function (e.g., with different finite terms, different partitioning of 
shared poles, different subleading or higher-order coiTections, or even a different ordering criterion). 

This, however, is not quite sufficient. Effectively, only the tree-level expansion of the shower 
would be affected by keeping track of such relative probabilities down along the shower chain; the 
Sudakov factors would remain unmodified. Such a procedure would therefore explicitly break the 
unitarity that is essential to resummation applications, leading to possibly exponentially different 
weights between the sets, which would be hard to interpret^. More intuitively, a big uncertainty 



'For example, two models that differ systematically by only a small amount on each branching, say 25%, would, after 
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on a veiy soft branching happening late in the shower should not be able to significantly change the 
entire event weight, jets and all. In the normal shower approach, it is the property of unitarity which 
keeps such things from happening; as soon as any correction grows large, its associated Sudakov 
factor must necessarily become small soon thereafter, keeping the total size of any correction inside a 
unit-probability integral. 

The main part of our proposal therefore concerns a simple way to restore unitarity explicitly also 
for the uncertainty variations, as follows. For each accepted branching, a number of trial branchings 
have usually first been generated and discai^ded, to eliminate the overcounting done by the trial func- 
tion. In ViNCiA, we have so far not been particularly careful to optimize the choice of trial function 
(see sec. 3.2), and hence we have quite many failed trials. These are relatively cheap to generate, how- 
ever, so the code is not significantly slowed by this inefficiency. Moreover, these failed trials actually 
turn out to be useful, even essential, in the present context. 

Just as eq. (123) expresses the relative probability for a branching to be accepted under two differ- 
ent sets of model parameters, 1 and 2, with 1 playing the role of phase-space generator and 2 the role 
of uncertainty variation, it is also possible to ask what the probability of a failed trial to have failed 
under different circumstances would have been. Thus each failed trial can actually be used to compute 
variations on the no-emission probability, i.e., the Sudakov factors. 

Specifically, for each trial, the no-emission probability for the model we use as our phase-space 
generator (which corresponds to the settings chosen by the user in ViNCiA, including matching, sub- 
leading corrections, etc.) is 

Pl;no = 1-Pl, (124) 

whereas the one for the alternative model should be 

P2;no = l-P2 = l-^Pl. (125) 

Osiai 

Thus, by multiplying the relative event weight W2/W1 by P2/P1 for each accepted branching and 
by -P2;no/-pL:no for cach failed one, we explicitly restore the unitarity of the set of weights {^2}- In 
order to prevent extreme outliers from substantially degrading the statistical precision of the variation 
samples, however, we limit the resulting weight adjustments to at most a factor of 2 per branching in 
the code (in either direction). The adjustment of the weights for the failed branchings takes the place 
of 'unfailing' those which should have succeeded with model 2. 

The accuracy of the approach obviously depends on the abundance of failed branchings. If the trial 
function is completely exact, and no branching ever fails, then the tree-level problem above will still 
occur. However, since ViNCiA typically generates significantly higher numbers of failed branchings 
than accepted ones, its effective numerical mapping of the changes in the Sudakov factors during the 
no-branching evolution periods should be reasonably accurate. 

To test whether the uncertainty bands produced in this way really reproduce what the shower 
model would have generated with different settings, we show a few distributions in Figs. 17 and 18, 
with default Vincia (thin blue line) plus an uncertainty variation (light blue band) on the left-hand 
side, and ViNCiA run with the actual settings corresponding to that variation on the right, for variations 
of the renormalization scale (Fig. 17) and of the antenna function finite terms (Fig. 18). In order to 
maximize the result of the variations, all matching is switched off, and hence the uncertainty bands 
are rather larger than would be the case for default ViNCiA settings. The L3 data (black points) [61] 
are included mostly to provide a constant reference across the plots; we postpone the discussion of the 



20 such branchings, differ by a factor 1.25^° = 100. If they differ by a factor of 2 instead, the result would be a million, 
clearly not a reasonable correction to the total event rate. 
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Figure 17: Comparison of unmatched results from VINCIA's automatic uncertainty variations in the 
Thrust observable around the default parameter set (left) with those from running the generator for 
each variation separately (right), for variation of the renomialization scale. The L3 data taken from 
ref. [61] is shown for comparison. The yellow band in the lower plots represents the experimental 
uncertainties on the thrust measurement. 

data comparison to sec. 8, where we discuss LEP observables. The top panels of each the plots shows 
MC compared to data, with both nonnalized to unity. The bottom panels show the ratio MC/data, with 
the uncertainties on the data shown as yellow shaded bands, the inner (lighter) one corresponding to 
the statistical component only and the outer (darker) shade corresponding to statistical plus systematic 
errors (added linearly, to be conservative). 

Comparing Figs. 17 and 18, one observes that the two different variations lead to qualitatively 
different shapes on the uncertainty predictions. The renormalization scale uncertainty. Fig. 17, pro- 
duces an uncertainty band of relatively constant size over the whole range of Thrust, whereas the finite 
terms. Fig. 18, only contribute to the uncertainty for large values of 1 — T, as expected. Comparing left 
to right in both figures, we conclude that both the features and the magnitude of the full uncertainty 
bands on the right are well reproduced by the weight variations on the left. 

Available Variations: So far, five types of automatic variations have been included in the Vinci A 
code, starting from version 1.025, via a simple on/off switch. These uncertainty variations are: 

• Vincia's default settings. This is obviously not a true uncertainty variation, but is provided as 
a useful comparison reference when the user has changed one or more parameters. 

• MAX and MIN variations of the renormalization scale. The default variation is by a factor of 2 
ai^ound p_\_. 

• MAX and MIN variations of the antenna function finite terms, as described in the online docu- 
mentation of the code^ . 

• Two variations in the ordering variable, one being closer to strong ordering in and the other 
to ordering in the m o variable. 

Available athttp://projects. hep forge .org/vincia/. 
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Figure 18: Comparison of unmatched results from VINCIA's automatic uncertainty variations in the 
Thrust observable around the default parameter set (left) with those from running the generator for 
each variation separately (right), for variation of the antenna-function finite terms. The L3 data taken 
from ref . [6 1] is shown for comparison. The yellow band in the lower plots represents the experimental 
uncertainties on the thrust measurement. 



• MAX and MIN variations of the subleading color corrections. The specific nature of the vari- 
ation depends on whether subleading corrections are switched on in the shower or not. If not, 
the MAX variation uses Ca for all gluon emission antennae and the MIN one C^?. If switched 
on, the coixection described in Section 4.3 is applied, but the correction itself is then modified 
by ±50% for the MAX and MIN variations here. 

These variations are provided as alternative weight sets for the generated events, which are available 
through methods described in the program's online manual. Limited user control over the variations 
is also included, such as the ability to change the factor of variation of the renormalization scale. 

When combining several variations to compute the total uncertainty, we advise to take just the 
largest bin-by-bin deviations (in either direction) as representing the uncertainty. We believe this is 
better than adding the individual terms together either linearly or quadratically, since the latter would 
have to be supplemented by a treatment of unknown conelations. With the maximal-deviation ap- 
proach, we are free to add as many uncertainty variations as we like, without the number of variations 
by itself leading to an inflation of the error. 

We should also note that, in the ViNCiA code, matching coefficients etc. are calculated for each 
uncertainty variation separately. The size of each band is therefore properly reduced, as expected, 
when switching on corrections that impact that particular- source of uncertainty. 

Finally, we note that, though the speed of the calculation is typically not significantly affected by 
adding uncertainty variations, the code does run slightly faster without them. We therefore advise to 
keep them switched off whenever they are not going to be used. 

7 Hadronization 

Since the ViNCiA code is a plug-in to Pythia, it is (almost) trivial to use Pythia's string hadroniza- 
tion model with ViNCiA, as long as one takes into account a few basic points: 
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• The matching of parton showers to hadronization models. 

• The "tuning" of the resulting shower+hadronization framework. 

Concerning the matching of parton showers to hadronization models, the main issue to keep in 
mind is that this matching is performed at a specific scale, the hadronization scale, which is imple- 
mented as a lower cutoff in the perturbative shower evolution, usually at a scale of order 1 GeV. Since 
no perturbative evolution is caixied out below that scale, the job of the hadronization model is then 
to give as good a representation as possible of all the physics that takes place at lower scales. Since 
the hadronization model is inherently non-perturbative, this means that the cutoff cannot be taken too 
high, or else it would become apparent that the hadronization model does not include a good descrip- 
tion of the perturbative parts. Vice versa, the cutoff cannot be taken too low, or else it would become 
apparent that the perturbative modeling does not include a good description of the non-perturbative 
parts. This is how one ends up with scales of order 1 GeV as the matching point. 

In principle, if the cutoff is varied around that point, both the perturbative shower and the non- 
perturbative modeling parameters should obey evolution equations that tell how each should scale, so 
that the end result would be approximatively independent of the cutoff. These evolution equations 
are nothing but an inclusive version of the shower evolution equations, which the shower obviously 
respects by definition. 

But there is so far no formalism for the non-perturbative modeling that allows us to take parameters 
"tuned" with one value of the cutoff and translate them for use with another value for the cutoff. Hence 
each setting of the parameters of the hadronization model are only valid for the exact cutoff value that 
they were tuned with. If one uses a lower cutoff, then there would be double counting between the 
shower, which now extends to lower scales, and the hadronization model, whose tuning attempted to 
absorb those same coixections as well as possible. Conversely, if one used a higher cutoff, there would 
be a kind of "dead zone" unreachable by evolution, and where also the hadronization model tuning 
did not attempt to absorb the corrections. 

To use Pythia's hadronization model directly with ViNCiA, we must therefore take the infrared 
cutoff to be at the same scale as the one used for the Pythia tuning, and since phase space is not 
one-dimensional, we also need to make sure it is in a variable which is as close to the one used by 
Pythia as possible. Contours corresponding to constant values of the Pythia 8 evolution variable 
were already illustrated in the discussion of evolution variables, fig. 8, and is reproduced in the left- 
hand pane of fig. 19 with an explicit cut showing the hadronization scale. 

In the dipole-antenna framework, it is not possible to work with exactly the same variable; since 
we do not keep the two sides of the dipole-antenna, a and h, separate, it is not possible to use a different 
form for the cutoff for the two sides, which is effectively what is done in Pythia. The closest we 
can come is to apply the cutoff if the smallest of the two scales, p^cvoi.a and p_Lcvoi,fe> is below the 
chosen cutoff scale, illustrated on the right pane of fig. 19. This will veto some branchings in ViNCiA 
that would have been allowed in Pythia, but since the Pythia radiation functions are not singular- 
in those regions, and since the kinematics of the corresponding phase-space points are near-colhnear, 
we do not expect this slight difference to have significant practical consequences. 

In addition to the possibility of using ViNClA's own variables, 2p^ or tud, as cutoffs, the option to 
use this emulation of the PYTHIA cutoff has therefore also been implemented in ViNCiA and should 
allow, to a first approximation, to use the Pythia hadronization model with any of the Vincia 
evolution settings, without retuning the non-perturbative parameters, as long as one accepts that the 
resulting answer will only be good up to perturbative uncertainties. That is, as long as the full ViNCiA 
uncertainty is estimated, the data should still be compatible with the resulting uncertainty bands. 
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Figure 19: Left: the hadronization cutoff in p±cvoi as it is imposed in Pythia, with the main plot 
showing one side of the evolving dipole and the inset the other. Right: in ViNCiA, the whole dipole- 
antenna evolves as one entity, and therefore the cutoff must be placed likewise. The remaining differ- 
ence only contains non-singular contributions from branchings in Pythia that accidentally throw the 
radiated paiton close to the recoiling one. 



We note that, in order to obtain a central ViNCiA tune, one would still have to perform a dedi- 
cated tuning of the Pythia 8 hadronization parameters using the particular Vinci A shower model 
for which the tuning is desired. This would also be necessary at the point when ViNClA's formal 
level of precision is higher than that of Pythia, in which case using hadronization settings tuned 
with Pythia's showers might actually result in incompatibility with the data, beyond the allowed 
perturbative uncertainty bands. 

A first step towards getting a dedicated ViNCiA tuning of the hadronization parameters was taken 
in April 2010 when one of the authors acted as host for a 1-week "industry internship" at CERN. 
Through the use of a specially developed runtime display for ViNCiA, and given some basic explana- 
tions about the effect of the different hadronization parameters on the LEP distributions, M. Jeppsson, 
a Danish middle-school student, succeeded in making a tune of ViNCiA to LEP data, including Thrust, 
the C and D parameters, jet rates, identified particle production rates, and the inclusive fractional mo- 
mentum distribution, 

2£'particlc 

a^particle = 1= ■ (126) 

The final parameters he settled on appeared well motivated and physical, and now constitute the 
default in ViNClA. They are given in appendix B for reference. The runtime display, which is based 
on ROOT, has subsequently been made publicly available as part of the ViNCiA package. To our 
knowledge, this study represents the first time "citizen science" has been used for event generator 
tuning. 



8 Comparison to LEP Data 

In the following, we have used version 1.025 of the ViNCiA plug-in and version 8. 145 of the Pythia 8 
generator, using default settings unless otherwise specified. Note that for ViNCiA, the default settings 
include a matching to tree-level matrix elements through third order in QCD (via its MadGraph 
interface), while Pythia only formally includes a first-order matching. 
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Figure 20: Comparison to the L3 light-flavor data set [61] (black points) at the Z pole for the 1 — T 
(left), C (middle), and D (right) event shape variables. ViNCiA is shown in thin blue lines, with 
shaded light-blue bands representing the perturbative uncertainty estimate. The middle pane on 
each plot illustrates the relative composition of the ViNCiA uncertainty band. For comparison, the 
PythiaS result is shown with a thick red line with open circles. The yellow bands in the bottom 
panels represent the experimental uncertainties on the measurement. 



To keep questions of mass effects separate (the implementation of which will be reported on in a 
separate paper [55]), we shall here mainly compare to a useful data set presented by the L3 collabora- 
tion [61], in which the contributions from light flavors (defined as u, d, s, c) has been sepai'ated from 
that of events containing b quarks. 

Unfortunately, however, the L3 light-flavor data set does not contain jet observables. We therefore 
include comparisons also to ALEPH and DELPHI jet observables that include all flavors, using a 
preliminary implementation of mass effects in ViNCiA [55]. Since the largest correction specific to b 
quarks is simply the B meson decay, for which we rely on Pythia's string hadronization and hadron 
decay model, we believe these comparisons are still meaningful, even if we must postpone a full 
discussion of them to the follow-up study in ref. [55]. 

In Fig. 20, we compare default ViNCiA and Pythia to the L3 light-flavor data for the Thrust 
(left) and the C (middle) and D (right) event shape parameters [61]. Dashed vertical lines indicate the 
boundaries between the 3- and 4-jet regions for the Thrust and C parameter (the right-most dashed 
line on the Thrust plot indicates the boundary with the 5-jet region). The D parameter measures the 
deviation from planar events and is a 4-jet observable over its entire range. Despite substantial dif- 
ferences in the shower modeling, matching level, and hadronization tune parameters, the two models 
give almost identical results. Further, since Pythia is already giving a very good description of this 
data, there is little for the additional matching in ViNCiA to improve on here. 

Still on Fig. 20, ViNClA's uncertainty bands give about a ±10% uncertainty over most of the 
observable ranges, with larger uncertainties near the edges of the distributions. The middle panels of 
the plots show the relative composition of the uncertainty estimates, and inform us that the renormal- 
ization scale variation is the dominant source of uncertainty for all the observables, with other sources 
only becoming competitive towards the right-hand extremes of the plots. This is an explicit conse- 
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Figure 21: Comparison to the L3 light-flavor data set [61] (black points) at the Z pole for the charged 
track multiplicity (left) and fractional momentum (right) spectra. ViNCiA is shown in thin blue lines, 
with shaded light-blue bands representing the perturbative uncertainty estimate. The middle pane on 
each plot illustrates the relative composition of the ViNCiA uncertainty band. For comparison, the 
PythiaS result is shown with a thick red line with open circles. The yellow bands in the bottom 
panels represent the experimental uncertainties on the measurement. 



quence of the tree-level matching in ViNCiA (by default imposed through third order), which signif- 
icantly reduces the allowed range of finite-term, ordering- variable, and subleading-color uncertainty. 
The renormalization-scale uncertainty, however, is unaffected by tree-level matching. Although the 
scale-dependence-reducing correction described in Section 4.1 is acting to reduce this dependence, 
the residual uncertainty from scale variation is still larger than that from any of the other sources. 

In Fig. 21, we compare to two infrared-sensitive observables also measured by L3 with light-flavor 
tagging, the charged track multiplicity (left) and the fractional momentum distribution (right), with 
the latter given by eq. (126). (Note that ref. [61] uses the notation = — Inx.) We conclude that 
Pythia 8 was probably tuned on slightly different observables, and hence the agreement obtained 
with ViNCiA is here improved both by giving a slightly naixower multiplicity distribution, with fewer 
low-multiplicity events and a slightly softer fragmentation spectrum, with fewer particles carrying x 
fractions very close to unity. One also notes that ViNClA's estimated uncertainty on the individual 
bins of the charged-track multiplicity distribution is much larger than the estimated uncertainty on the 
fragmentation spectmm. Recall, however, that ViNCiA is only able to vary the perturbative parameters 
— variations of the string fragmentation parameters would have to be included here to gain a better 
understanding of the full uncertainties. All we can say at this level is that the charged-multiplicity 
distribution appears to suffer from a larger perturbative uncertainty than the fragmentation spectrum. 

A further set of variables that is interesting in the context of differential multi-jet production 
are the so-called four-jet angles, which were also measured at LEP. Not having found a public data 
repository containing this particular- data, however, we instead resorted to extracting the data point 
values from the Herwig-i-i- source code [39], where it is encoded for validation and tuning purposes. 
A comparison between this data and default ViNCiA and PYTHIA is shown in Fig. 22. Again, it is 
clear that Pythia itself is already doing a very good job. Since Pythia is not matched to 4-jet 
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Figure 22: Comparison to DELPHI 4-jet angle measurements (black points) at the Z pole. ViNCiA 
is shown in thin blue lines, with shaded light-blue bands representing the perturbative uncertainty 
estimate. The middle pane on each plot illustrates the relative composition of the ViNCiA uncertainty 
band. For comparison, the PythiaS result is shown with a thick red line with open circles. The 
yellow bands in the bottom panels represent the experimental uncertainties on the measurement. 



matrix elements and also does not contain explicit spin correlations in the shower, this may at first be 
surprising. However, Pythia does coiTclate the production and decay planes of gluons in the shower, 
and thereby includes the leading effect of gluon polarization. The ViNCiA shower, on the other hand, 
contains no polarization effects a priori. In ViNClA's case, the effective correlations of the four-jet 
angles are instead coming from matching to the 4-parton matrix elements, and both codes are able 
to describe the 4-jet angles within a roughly 5% margin, which is comparable to the experimental 
precision. 

Finally, in Fig. 23, we compare to the jet resolutions measured by the ALEPH experiment [62]. 
Firstly, note that pure Pythia is basically able to describe all the distributions, within the experimen- 
tal accuracy, despite its being matched only to Z — 3 partons. On the one hand, this is good, since it 
implies that the Pythia 8 shower is delivering a quite good approximation to QCD also beyond the 
matched orders. On the other hand, it also means that we are not really able to quantify any signifi- 
cant improvement by matching using this data alone. This may partly be due to the data having quite 
large uncertainties for hard multi-jet configurations; at least ±20% for jet resolutions ln(y45) > —4, 
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Figure 23: Comparison to ALEPH jet resolution measurements [62] (black points) at the Z pole. VlN- 
CIA is shown in thin blue lines, with shaded light-blue bands representing the perturbative uncertainty 
estimate. The middle pane on each plot illustrates the relative composition of the ViNCiA uncertainty 
band. For comparison, the PythiaS result is shown with a thick red line with open circles. The 
yellow bands in the bottom panels represent the experimental uncertainties on the measurement. 
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which coiTesponds to p± scales of order 10 GeV. Thus, this data does not really have enough power to 
make precision tests in the hard multi-jet region. This relatively low discriminating power of the hard 
multi-jet data will likely be an even greater problem for testing any future one-loop matched schemes. 

We see four possible avenues to improve on this situation, requiring various levels of additional 
experimental effort. One, combine the hard multi-jet data from the four LEP experiments to increase 
the overall discriminating power. Two, re-process the data focussing on observables more specifically 
tailored to project out regions not dominated by leading logs, e.g., by measuring Unn+i while impos- 
ing specific constraints on the resolutions of the other (n — 1) hai^der jets in the event. Thr^ee, increased 
statistics from a future facility such as GigaZ. Four, measurements of hard jet rates and resolutions 
outside the e^e~ environment, i.e., at hadron colliders, again possibly placing more exclusive restric- 
tions on the multi-jet stmcture of the events. Especially the latter looks promising now in the dawn of 
the LHC era, but of course comes at the price of introducing additional uncertainties due to the colored 
initial states. LEP is therefore likely to retain its position as our main jet fragmentation laboratory for 
the foreseeable future, and with the official closing of the LEP experiments this year, we wish to en- 
courage those in a position to do so to keep the LEP data 'alive and well' for future analysis studies 
that are likely to involve tests of models far more sophisticated than those that were available ten or 
twenty years ago. This makes it necessary to look at data much more differentially and/or exclusively 
than is possible with, e.g., with the observables that were included in our comparisons here. 

9 Conclusions 

We have taken the next step in developing the formalism, started in ref. [33], for generic parton 
showers based on the dipole-antenna formalism. Evolution equations for a wide class of evolution 
variables, kinematics maps, and radiation functions, have been presented, including all the necessary 
steps to construct an explicit stochastic Markov-chain Monte Carlo code. 

The basic ideas behind this shower model are similar to those behind the existing Ariadne pro- 
gram [17], to whose properties we make some comparisons. Aside from the more generic formalism, 
we also propose some systematic improvements, including suppressed unordered branchings to cover 
the hard region of phase space and systematic "next-to-leading-color" (NLC) corrections. We com- 
pare explicitly to matrix elements at both leading and subleading color for Z — 4, 5, and 6 partons, 
to check the validity of our approximations. 

We have also presented a new method for matching parton showers to tree-level matrix elements 
at the multi-jet level, formulated in a language appropriate to our shower framework. At lowest 
order, it is similar to an older scheme by Sjostrand and collaborators [28, 29], which has recently 
been reformulated in a more generic NLO context called POWHEG [30]. Though our scheme is 
therefore similar to, and compatible with, these existing methods, we here extend the method to tree- 
level matching involving more than one emission. As such, the method is at the same formal level 
of precision as the Menlops approach [26], but with the difference that the multi -parton matrix 
element coiTections are here exponentiated, which should both improve the logarithmic accuracy of 
our Sudakov factors at the same time as making the approach much less sensitive to so-called matching 
scales (scales below which the matrix-element corrections are switched off). We do still advocate 
imposing some matching scale at high multiplicities, mostly to avoid spending lots of time computing 
matrix elements for very soft emissions that the unmatched shower is describing coiTectly anyway. 
For the results reported on in this paper, a matching scale of 5 GeV (above the hadronization scale but 
well below typical jet pj'S) was therefore imposed starting from third order in QCD (corresponding to 
Z — )• 5 partons). 
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By default, the generated events form one continuous sample with no specific separation between 
jet multiplicities and with all events having unit weight. 

The matched shower algorithm is implemented in the ViNCiA plug-in to the Pythia 8 Monte 
Carlo generator. Apart from its dependence on Pythia 8, the plug-in is self-contained with its own 
documentation and user-definable parameters. It also includes facilities for linking to the FastJet 
[63] and ROOT packages. In the latter case, a mn-time interface has been developed that allows to 
display ROOT histograms in real time during the generator run, which can be useful both to give an 
immediate sanity check that histograms are being filled correctly, and also to visualize the gradual 
improvement in MC statistics over the run. The generated events have similar physical properties 
as those generated by standai^d Pythia 8 and can be passed through the latter's string model for 
hadronization, and, e.g., to HEPMC [64] for further processing, e.g., by analysis tools like RiVET [65] 
or detector simulation packages. 

Finally, we have presented a new efficient and automatable method for the evaluation of uncer- 
tainty bands by parton shower generators. The method draws on the unitarity property of shower 
calculations to compute several sets of weights for a single generated event sample. This sample can 
then be subjected to cuts, hadronized, passed through detector simulations, etc., and the uncertainty 
variations can be obtained by filUng histograms with each of the different weight sets separately at 
any stage during the processing. We have implemented this method in ViNCiA, to provide an option 
for automatic evaluation of renormalization-scale, finite-tem, ordering, and subleading-color uncer- 
tainties. 
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A Comparisons to 2nd Order QCD Matrix Elements 



In the following, we present a series of figures comparing dipole-antenna showers — ordered in 
various different evolution variables — to matrix elements, on the R4 ratio defined in eq. (83). As in 
sec. 3, these plots were made on 20M random points in a flat phase-space scan, using Rambo [50]. 
We use the default ViNCiA antenna functions and kinematics maps in all cases and only vaiy the 
evolution variable. 

On each plot, the quantity on the x axis is a scale characterizing the first emission (2 — 3), while 
the quantity on the y axis is a scale characterizing the second (3 — )• 4). Further, the quantity on the 
y axis has been normalized to the one on the x axis, so that it really represents the ordering of the 
second emission, relative to the first. Thus, phase-space points with hard initial 2 — 3 branchings 
will lie to the right in the figures, while soft (strongly ordered) initial emissions lie to the left. On the 
y axis, points where the second scale is much smaller than the first (i.e., strongly ordered relative to 
the first) will be towards the bottom of the plot, while unordered points will be towards the top of the 
plot. 

To help the eye, we have added a horizontal red line at ln(y) = on the plots, dividing the phase 
space for the 2nd emission into an ordered part (below the line) and an unordered one. Similarly, the 
emphasized black box highlights the region where x is more than an order of magnitude below Mz, 
and y more than an order of magnitude below x, i.e., the doubly-ordered region. Finally, the dotted 
Unes show contours of constant ln(y/M|), or equivalently constant y. 

For each evolution variable, we show 4 plots, all with logarithmic x and y scales: 

• Top Left: Average PS/ME ratio {R4) vs. the p± scales of the two emissions. On the y axis, the 
smallest of the two possible pj_ scales in the 4-parton configuration. On the x axis, the p± scale 
of the corresponding 3-parton configuration. 

• Top Right: The RMS width of the left-hand ratio, counting dead zones as having a factor 10 
difference (otherwise the log of the ratio would be undefined). This helps illustrate whether a 
good average agreement on the left-hand side is just an accident of a wide distribution centered 
around unity, or whether the distribution itself is really naiTow. 

• Bottom Left: Same as top left, but with the y axis showing the invariant mass of the two gluons. 
This projection is interesting since it allows us to isolate a particular double-collinear limit, as 
follows. If rn^g ~ m^, then the two gluons are well separated and carry all the energy of the 
original paitons. This energy can only have been transfeiTcd to them by two successive extreme 
coUinear branchings, one where the quark gives all its energy to the first gluon, and a second 
where the antiquark gives all its energy to the other gluon. To check if we get this non-trivial 
limit right, the upper edge of this distribution is therefore especially interesting. 

• Bottom Right: Same as top left, but with the ordering measure being invariant mass instead of 
p±. This distribution is complementary to the top left one, showing the same ratio in a slightly 
different projection, versus a measure of invariant masses, again in order to check whether any 
agreement observed in the above variables persists when doing a different cut thi^ough phase 
space. 

Note that, since we use the default Vincia antenna functions, which reproduce the Z — )• 3 matrix 
elements, we expect a good agreement even when the first emission is not strongly ordered. It is 
therefore mainly the ordering of the second emission with respect to the first that is interesting. 
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Phase-Space Ordering 
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Figure 24: Phase-space-ordered antenna approximation compared to 2nd order QCD matrix elements. 
Note: this roughly corresponds to a mass-ordered parton shower without coherence. Although the 
double-soft limit is eventually reached, there is a large overcounting over most of phase space, reflect- 
ing a lack of coherence. Also, the double counting extends into the double-collinear region at the top 
of the lower left-hand plot. This ordering, therefore, does not lead to the correct multiple-coUinear 
singular limit. 
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Transverse-Momentum Ordering (ARIADNE) 
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Figure 25: Transverse-Momentum-Ordered antenna approximation compared to 2nd order QCD ma- 
trix elements, using the Ariadne definition of p±, wiiicii is also the default evolution variable in 
ViNClA. Most of the double-counting evident for phase-space ordering has been removed, and the 
shower approximation now also gives the correct answer in the double-collineai^ region at the top of 
the lower left-hand plot. The price is the introduction of a dead zone, visible at the top of the upper 
left-hand plot. The size of the dead zone in the flat phase-space scan amounts to about 2% of all 
sampled points. 
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Transverse-Momentum Ordering (PYTHIA) 
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Figure 26: Transverse-Momentum-Ordered antenna approximation compared to 2nd order QCD ma- 
trix elements, using Pythia's definition of Transverse Momentum, p^cvoi- Note: the antenna func- 
tions and kinematics maps are still the default ViNCiA ones, hence these results do not directly corre- 
spond to what would be obtained with the PYTHIA program. The PYTHIA p^-definition is a bit more 
restrictive and increases the size of the dead zone to 5% of the sampled points. This in turn leads 
to a lower average ratio in some regions. The fact that the RMS distribution indicates a large spread 
of weights even on the edges of the strongly-ordered region (black box) reflects the migration of a 
few dead-zone points to this region, due to the mismatch between the ordering variable and the p± 
definition on the axes. 
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Mass-Ordering (VINCIA) 




Figure 27: Mass-Ordered antenna approximation compared to 2nd order QCD matrix elements, using 
ViNClA's definition of Daughter Dipole Mass, m/j. This ordering is a bit less restrictive than p±, 
hence the size of the dead zone shrinks to about 1 % of the sampled points. The tradeoff is that some 
overcounting remains over some parts of phase space. As for the Pythia -ordering, the RMS 
distribution indicates a quite wide distribution extending inside the doubly-ordered box, which we 
interpret as being due to the mapping of some dead-zone points to this region. 
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Energy-Ordering (DM) 
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Figure 28: Energy-Ordered antenna approximation compared to 2nd order QCD matrix elements, 
using a definition of energy a la Doksliitzer-Marchesini (DM). Although a small dead zone in the 
unordered region still exists (0.6% of the sampled points), there remains a very large overcounting 
over significant parts of phase space, including the double-collinear region mentioned before, at the 
top of the lower left-hand plot. We conclude that this variable does not lead to the correct multiple- 
coUinear singular- limit. 
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Modified Energy-Ordering (VINCIA) 
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Figure 29: Energy-Ordered antenna approximation compared to 2nd order QCD matrix elements, 
using a modified definition of energy that makes it explicitly sensitive to collinear emissions. The 
dead zone is still quite small, ca. 0.7% of the sampled points, but the approximation is improved with 
respect to energy ordering over most of phase space, including the double-collinear region at the top 
of the lower left-hand plot. 
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Angular-Ordering (HERWIG++) 
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Figure 30: Angular-Ordered antenna approximation compared to 2nd order QCD matrix elements, 
using the HERWIG++ definition of angles. Note: the antenna functions and kinematics maps are still 
the default ViNCiA ones, hence these results do not directly con^espond to what would be obtained 
with the HERWIG++ program. The dead zone in this variable amounts to 1% of the sampled points. 
Although the double-coUinear region appears to be well approximated (see top of lower left-hand 
plot), there appears to be a tendency toward untercounting of the double-soft region (box in top left 
plot) . The RMS spread is also very lai^ge over substantial regions of phase space (top right plot), 
extending well beyond the region affected by the dead zone. 
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V-Ordering (VINCIA) 




Figure 31: Angular-Ordered antenna approximation compared to 2nd order QCD matrix elements, 
using Vincia's modified angular ordering, Vs, which is explicitly sensitive to soft and coUinear emis- 
sions. The overall picture is quite similar to that for the Herwig-i-i- definition of angular-ordering. 
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Smooth Transverse-Momentum Ordering (VINCIA) 
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Figure 32: Transverse-momentum-ordered antenna approximation compared to 2nd order QCD ma- 
trix elements, using Ariadne's definition of pj_ and Vincia's smootfi suppression factor instead of 
the usual strong ordering condition. This corresponds to the default in ViNCiA without matching. 
(Note: by default, matching to Z — )• 4 is on in ViNCiA, over all of phase space, and hence these ratios 
are all equal unity). 
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Smooth Mass-Ordering (VINCIA) 
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Figure 33: Mass-ordered antenna approximation compared to 2nd order QCD matrix elements, us- 
ing Vincia's definition of Md applied as a smooth suppression factor instead of the usual strong 
ordering condition. Although the dead zone has been removed without introducing a catastrophic 
over-counting, the results are less impressive than for p^-ordering. Note: the default in ViNCiA is 
therefore to apply the smooth suppression factor in p± regardless of the actual choice of evolution 
variable. 
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B Jeppsson Tune Parameters 



The following table gives an overview of the parameters used for ViNCiA and Pythia for the results 
obtained in this paper. The default settings in Pythia 8.145 were obtained by a fit to a large amount 
of LEP data, using the Rivet+Professor framework [65,66]. The Vincia 1.025 parameters were 
tuned manually, as reported on in sec. 7. The latter included the total charged particle multiplicity and 
^ = — log(x) distributions, mean multiplicities of light mesons and baryons, event shapes (1 — T, 
Major, Minor, C, D, Oblateness), and jet resolution scales (y23^ 2/34^ 2/45^ yse) extracted from the 
measurements contained in the HepData [67] repository. 



Parameter 


Pythia 8.145 Default 


Vincia 1.025 


Shower Parameters 


Evolution Variable 


P±evol 


Pi. 


Renormalization Scale 


P±evol 


P± 


as{Mz) 


0.1383 


0.139 


as loop order 


1 


1 


Shower Cutoff Variable 


P±evol 


2p± 


Shower Cutoff Scale in GeV 


0.5 


1.0 


String Breakup Parameters 


StringZ : aLund 


0.30 


0.28 


StringZ : bLund 


0.80 


0.55 


StringPT : sigma 


0.304 


0.275 


StringPT : enhancedFraction 


0.01 


0.01 


StringPT : enhancedWidth 


2.0 


2.0 


Flavor Parameters 


StringFlav: probStoUD 


0.19 


0.21 


StringFlav: mesonUDvector 


0.62 


0.40 


StringFlav: me sons vector 


0.725 


0.6 


StringFlav: probQQtoQ 


0.09 


0.079 


StringFlav :probSQtoQQ 


1.0 


1.0 


StringFlav :probQQltoQQO 


0.027 


0.035 


StringFlav: decupletSup 


1.0 


1.0 


StringFlav: etaSup 


0.63 


0.60 


StringFlav: etaPrimeSup 


0.12 


0.075 
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